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1.0  INTRODUCTION 


Effective  transonic  operation  is  required  of  modern-day  transport  and 
fighter  aircraft.  This  consideration  has  stimulated  development  of  numerical 
methods  to  compute  transonic  flowfields  about  increasingly  complex 
geometries.  In  order  to  be  useful  for  design  and  performance  analysis,  such 
methods  must  contain  accurate,  reliable  schemes  for  solving  the  nonlinear 
equations  governing  the  flowfields  as  well  as  the  capability  to  handle 
realistic  configurations. 

Various  finite-difference  methods  have  been  developed  to  calculate  three- 
dimensional,  transonic  potential  flowfields  about  isolated  wings  and  wlng/body 
configurations.  To  date,  these  methods  have  used  formulations  that  require 
alignment  of  rjmputlng  coordinates  with  appropriate  geometric  surfaces  in 
order  to  apply  surface  boundary  conditions.  The  codes  of  Ballhaus  and 
Bailey^”^  and  Boppe,'^”^  which  solve  various  forms  of  the  transonic  small- 
disturbance  potential  equation,  use  mean-plane  and  nearest-point 
approximations  to  represent  wing  and  fuselage  geometries,  respectively;  the 
approximating  surface  is  chosen  to  coincide  with  a  cartesian  grid  and  is  the 
location  where  a  linearized  f low-tangency  boundary  conditions  is  enforced. 

The  codes  of  Jameson  and  Caughey^"^  are  based  on  the  full  potential  equation 
and  require  application  of  an  exact  surface  condition  using  analytic  and 
numerical  mappings  which  generate  computational  grids  that  conform  to  the 
entire  surface  geometry.  Although  the  small-disturbance  approach  has  so  far 
produced  the  most  extensive  geometry-handling  capability , comparative 
calculations  indicate  that  methods  based  on  the  full  potential  equation 
provide  greater  accuracy  in  matching  experimental  data.^® 

This  report  presents  an  alternative  approach  to  transonic  wing/body 
calculations  in  which  the  full  potential  equation  is  solved  using  coordinates 
that,  in  general,  do  not  coincide  with  the  configuration  surface.  The 
wing/body  surface  appears  as  a  curved  boundary  within  a  rectangular  computing 
grid.  In  order  to  enforce  the  f low-tangency  condition  exactly  on  this  curved 
boundary,  an  Imaging  scheme  is  used  which  Involves  surface  control  points  and 
surface-adjacent  grid  points.  At  the  expense  of  some  computational  detail, 
the  approximations  inherent  in  a  small-disturbance  formulation  and  the 
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complicated  mappings  required  to  produce  a  surface-conforming  grid  are  both 
avoided.  Schemes  similar  to  the  one  outlined  have  been  applied  previously  to 
calculate  flows  about  simpler,  two-dimensional  planar  and  axisymmetric 
conf igurations , ^^-14 

The  objective  of  this  contract  was  to  develop  a  computer  program,  based 
on  a  global  system  of  non-surface-fitted  coordinates,  to  calculate  transonic 
flowfields  about  wlng/f Inite-fuselage  configurations  and  to  compare  calculated 
solutions  with  relevant  experimental  data.  It  was  determined  part  way  into 
the  contract  period  that  the  coordinate  formulation  originally  proposed  as  the 
basis  for  this  effort  was  Inadequate  and  a  more  complicated  system.  Involving 
a  higher  degree  of  nonorthogonality,  was  needed  to  achieve  acceptable 
solutions.  The  delay  associated  with  this  reformulation  as  well  as  some 
initial  difficulty  in  achieving  stable  operation  of  the  program  have  resulted 
in  a  geometry  capability  that  consists  of  a  general  wing  attached  to  fuselage 
of  seml-lnf Inlte  length. 

Sections  2.0  and  3.0  of  this  report  outline  the  formulation  of  the 
wing/body  problem.  Section  4.0  describes  the  numerical  solution  scheme.  A 
discussion  of  results  is  given  in  Section  5.0,  and  remarks  relating  to 
required  code  improvements  and  possible  future  directions  are  contained  in 
Section  6.0.  Appendix  A  is  a  summary  user's  guide  for  the  computer  program, 
and  Appendix  B  is  a  program  listing. 
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2.0  GEOMETRY  AND  COORDINATES 


Consider  the  configuration  shown  in  Figure  1.  A  wing  with  raultiple 
sweep,  taper,  and  dihedral  is  attached  at  a  mid-height  position  to  a  blunt- 
nosed  cylinder  of  semi-infinite  length.  The  wing  leading  and  trailing  edges 
are  comprised  of  straight-line  segments;  each  edge  may  (but  need  not)  contain 
a  single  break  at  a  spanwise  location  that  differs  from  the  other.  The  wing 
varies  arbitrarily  along  its  span  in  section  profile  and  section  twist 
angle.  It  is  assumed,  however,  that  the  variation  of  wing  properties  is 
piecewise  linear  between  defining  sections.  The  fuselage  component  of  the 
configuration  has  an  asymmetric  side  profile  and  a  circular  crossplane  section 
of  varying  radius  along  its  length. 

For  convenience  in  formulating  the  numerical  scheme  to  calculate  the 
flowfleld  about  this  class  of  configurations,  a  global  system  of  wing-adapted 
coordinates  Is  introduced  followed  by  a  series  of  coordinate  stretchings. 

This  accomplishes  several  purposes.  The  wing/body  geometry  is  transformed  to 
a  standard  and  relatively  simple  form  to  ease  computational  bookkeeping. 

Also,  the  infinite  physical  domain  is  mapped  to  a  finite,  rectangular 


paralleloplped  along  whose  outer  boundaries  exact  far-field  boundary 
conditions  can  be  applied*  (However,  the  transformed  wing/body  surface 
remains  embedded  within  the  parallelepiped  as  an  irregular,  interior 
boundary.)  Finally,  grid-point  spacing  can  be  arranged  to  provide 
satisfactory  surface  resolution,  at  least  on  the  wing  at  the  present  stage  of 
development,  as  well  as  efficient  distribution  of  points  throughout  the 
computational  domain  (dense  near  the  configuration  surface  and  progressively 
more  sparse  as  distance  from  the  surface  Increases). 

The  wing-adapted  coordinates  are  defined  by  the  following  relations; 


X  -  X,  (y) 


Y  =  Z 
^  b  • 


z  -  z^(x,y) 
T(y)  c(y) 


»here  X2^g(y)  represents  the  streamwise  position  (sweep)  of  the  wing  leading 
edge,  c(y)  is  the  wing-section  chord,  b  is  the  wing  span,  r(y)  is  the  wing- 
section  thickness-to-chord  ratio  normalized  by  its  value  at  the  wing  root,  and 
Zj^(x,y)  is  a  nonplanar  mean  surface  determined  by  the  wing  geometry.  This 
mean  surface  is  specified  as  follows; 


\(x,y)  «  z^^(y) 

[x  < 

Xig(y)], 

(2a) 

=  ^ig(y) 

-  [x  -  x^^(y)]  tan  ot^(y) 

[Xie(y)  <  X  < 

Xte(y)K 

(2b) 

=  ^ie(y) 

-  c(y)  tan  a^(y) 

[x  > 

Xtg(y)]. 

(2c) 

where  denotes  the  vertical  position  (dihedral)  of  the  wing  leading 

edge,  oij.(y)  is  the  wing-section  twist  angle,  and  Xj.g(y)  gives  the  streamwise 
location  of  the  wing  trailing  edge.  Note  that  c(y)  =  -  x^g(y). 

The  transformation  Equations  (1)  effectively  shear  out  wing  sweep,  taper 
in  both  chord  and  thickness,  dihedral,  and  twist.  In  the  (X,Y,Z)  coordinates, 
the  wing  appears  to  be  planar  with  a  rectangular  planform  and  uniform 
thickness  over  its  entire  span.  The  constant  in  the  x-transformation  shifts 
the  X-origin  to  the  wing  mid-chord  locus;  the  leading  and  trailing  edges  are 
at  X  »  ±  1/2,  respectively.  Lines  of  constant  X  coincide  with  constant- 
percentage-chord  stations  along  the  span.  The  Y  origin  is  at  the  wing/body 
centerline,  and  the  wingtip  is  at  Y  =  1/2.  The  Z-origln  is  on  the  mean 
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surface  Zj^(x,y).  Except  for  the  nonplanar  wing  mean  surface  the 

(X,Y,Z)  coordinates  are  similar  to  those  used  in  transonic  small-disturbance 
formulae  ions . 

The  spanwise  region  containing  the  fuselage  and  the  region  outboard  of 
the  wingtip  are  treated  by  extending  the  wing  leading  and  trailing  edges  while 
holding  the  section  twist  angle  fixed  at  its  wing  root  and  wingtip  values, 
respectively.  The  fuselage  transforms  in  a  regular  but  not  well-defined 
manner  with  distortion  in  length,  profile,  and  crossplane  section.  In  dealing 
with  the  outboard  region,  the  normalizing  chord  length  in  Equation  (1)  is  held 
fixed  at  a  constant  value  efy^)  beyond  an  arbitrarily  specified  spanwise 
position  y  =  y2  to  prevent  crossing  of  the  extended  wing  edges. 

In  order  to  produce  a  finite  computational  domain,  each  of  the  (X,Y,Z) 
coordinates  is  stretched  individually  as  described  below.  The  stretching 
formulas  are  adaptations  to  the  present  application  of  functions  used 
successfully  to  compute  the  transonic  flow  about  an  airfoil.  Figure  2  shows 
the  relationship  between  the  chordwise  and  spanwise  coordinates  and  also 
indicates  the  different  regions  of  the  physical,  intermediate,  and  stretched 
domains . 


•  X-coordinate  to  f-coordlnate 


Region  II:  X 

Region  1,111;  X 


f  (aj  +  32 

T  X  +  A  tan  r^(f.  ±  i  ) 
o  1  2  o 


+  A-  tan 


(3a) 

{3b) 


In  Equation  (3b),  the  upper  set  of  signs  refers  to  region  1  and  the  lower  set 
to  region  III.  The  three  streamwise  regions  encompass  the  following  ranges: 


1:  ^  <  X  <  -X  ,  -(1  +  t  )  <  T  <  -r  , 

o  o  o 


II:  -X  X  <  X  , 
o  o 


[II:  X  <  X  <  , 

o 


-r  <  f  <  r 
o  O 


<-  <  <  (1  +  ':  ) 
o  o 


(4a) 

(4b) 

(4c) 


.>0 

t 

Ob. 

i 

K  (y ) 
le 

00 

X  0 

‘  00 

» 

0 

0 

►  ^  ) 

t  0 

( 1  •  t  1 
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r* 
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(a)  Wing-section  plane 

I  -  I  I  -I  - -  rj 
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(b)  Wina/body  crossplane 
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Figure  2.  Coordinate  reialionsbips. 


Regions  1  and  III  nominally  represent  the  regions  upstream  of  the  wing  leading 

edge  and  downstream  of  the  trailing  edge,  respectively,  while  Region  II  covers 

the  wing-section  chord.  This  stretching  is  symmetric  about  the  origin. 

Constants  aj^  and  a2  are  determined  by  the  conditions  X  =  and 

dX/d£  =  •irA,/2  at  C  =  f  .  The  transition  stations  ±X  between  the  cubic  and 
1  o  o 

tangent  stretching  functions  are  chosen  to  occur  a  small  distance  inside  the 
wing  leading  and  trailing  edges.  Parameter  determines  how  much  of 
the  f -domain  is  confined  between  the  wing  edges.  For  an  evenly  spaced 
f-grld,  the  number  of  £-steps  in  Region  II  compared  with  the  total  number 
of  f-steps  will  be  in  the  ratio  C^^/d  +  C^).  Constants  Aj  and  A2  mainly  allow 
control  over  grid-point  spacing  in  regions  I  and  III  but  also  provide  some 
control  over  the  uniformity  of  spacing  in  Region  II. 
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•  Y-coordinate  to  n-coordinate 


Region  IV;  Y  =  n  (b  (5a) 

Region  V:  Y  -  Y^  +  tan  ~  ^  ^2  ~ 

Spanwise  regions  cover  the  following  ranges: 

IV:  0  <  Y  <  Y  ,  0  <  n  <  n  ,  (6a) 

o  o 

V:  Y<Y<“,  n<n<(l+n).  (6b) 

o  o  o 

This  spanwise  stretching  is  similar  to  the  one  used  in  the  streamwise 
direction  but  is  applied  only  to  the  half-space  because  only  unyawed  wing/body 
configurations  are  considered.  Region  IV  extends  over  the  wing/body  semispan, 
and  Region  V  encompasses  the  domain  outboard  of  the  wingtip.  Constants  b^  and 
b2  follow  from  the  requirements  that  T  =  Y^  and  dY/dn  =  ■nBj/2  at  n  =  h^.  The 
transition  station  Y^  is  set  slightly  inside  the  wingtip,  and  parameter 
determines  the  fraction  of  the  n-domain  confined  to  the  semispan  region. 
Constants  Bj  and  Bo  control  spacing  of  the  grid  outboard  of  the  wing  tip,  in 
region  V. 

•  Z-coordinate  to  ^-coordinate 

Z  =  C^  tan  (iiC/2)  (7) 

Constant  Cj  controls  vertical  grid-point  spacing  near  the  wing  mean  surface. 
The  range  of  this  stretching  is: 

-«  <  Z  <  “>  ,  -1  <  C  <  1  .  (8) 

Figure  3  is  a  schematic  representation  of  the  nonorthogonal  (^,n,C) 
coordinate  system  which  indicates  its  relationship  to  the  wing  geometry. 

Figure  4  shows  a  schematic  representation  in  the  physical  (x,y,z)  domain 
of  a  grid-point  distribution  that  is  uniformly  spaced  in  each  of  the  (C,h,4) 


directions.  The  configuration  shovm  involves  a  planar  wing  with  a  straight 
leading  edge.  Each  point  in  the  schematic  actually  represents  a  line  of  grid 
points.  Grid  taper  and  the  clustering  of  points  near  the  wing  mean  surface 
and  between  the  wing  edges  are  clearly  evident.  The  dashed  line  corresponds 
to  the  spanwise  station  y  =  discussed  previously.  The  planform  and 
crossplane  views  of  Figure  4  emphasize  the  wing  orientation  of  the  present 
coordinates  and  illustrate  their  lack  of  suitability  for  representing  the 
fuselage  region.  A  method  for  improving  resolution  about  the  fuselage  is 
described  in  Section  5.0.  The  wing-section  detail  emphasizes  the 
nonconformity  between  the  computation  coordinates  and  the  geometry  surface. 

For  later  use  in  writing  the  transformed  differential  equations  that 
govern  the  wing/body  flowfield,  the  stretching  functions  are  denoted 
symbolically  as 

4  =  4(X),  n  =  n(Y),  ;  =  4(Z).  (9) 

Stretching  derivatives  are  then  defined  bv  the  relations 

f(4)  =  d4/dX,  g(n)  =  dn/dY,  h(;)  »  dC/dZ.  (10) 

Evaluation  of  these  derivatives  follows  from  Equations  (3),  (5),  and  (7), 
respectively. 
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3.0  GOVERNING  EQUATIONS 


3 • I  Physical  Domain 

The  steady,  inviscid,  transonic  flow  past  a  wing/body  configuration  can 
be  considered  isentropic,  and  hence  irrotational,  if  only  weak  shock  waves  are 
present.  Such  a  flow  can  be  characterized  by  a  velocity  potential  $(x,y,z) 
that  is  related  to  the  streamwlse  x-component,  the  spanwlse  y-component,  and 
the  vertical  z-component  of  velocity  by  the  relations 


u='I>,v=<li,w=$,  (11) 

X  y  z 

where  coordinate-symbol  subscripts  denote  partial  differentiation.  In 
Equation  (11)  and  below,  it  is  assumed  that  all  velocities  are  normalized  by 
the  freestream  speed  . 

Since  $  is  singular  at  infinity,  it  is  convenient  to  Introduce  a 
perturbation  potential  ((i(x,y,z)  according  to  the  expression 

4>  =  X  cos  a  +  z  sin  a  +  4)  ,  (12) 

where  a  is  the  angle  of  attack  of  the  incident  flow  measured  in  a  wing-section 
(x-z)  plane.  Then,  the  governing  equation  for  the  flowfield,  written  in 
cartesian  coordinates,  has  the  form 

(a^  -  u^)(t)  +  (a^  -  v^)^  +  (a^  ~  w^)4i  -  2uv4> 

XX  yy  xy 

-  2vw4i  -  2uw4i  =0,  (13) 

yz  xz 

where 


u  “  cos  a  +  (p  , 

X 


(14a) 


V  =  , 


(14b) 


w  “  sin  a  +  4> 


(14c) 
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The  local  speed  of  sound  a  is  determined  by  the  relation 


00 


0  ')  ')  ') 

where  q'^  =  u  +  v  +  w  ,  is  the  freestream  Mach  number,  and  y  represents 

the  ratio  of  specific  heats  of  the  medium  (for  air,  y  =  1.4). 

The  pressure  coefficient  at  any  point  in  the  flowfield  is  given  by  the 
expression 


C 

P 


Y/(Y-1) 


(16) 


The  boundary  condition  for  Equation  (13)  on  the  wlng/body  requires  that 
the  flow  be  tangent  to  the  surface  and  can  be  written  as 


(cos  a  + 


//'  +  -  (sin  a  + 

X  y  y 


0; 


surface 


(17) 


derivatives  and  are  streaniwise  and  spanwise  surface  slopes, 
respectively.  In  addition,  the  Kutta  condition  requires  that  a  circulation 
r(y^)  must  exist  at  each  spanwise  wing  station  y^  which  is  of  such  magnitude 
that  the  flow  passes  smoothly  off  the  sharp  trailing  edge.  A  vortex  sheet 
extends  behind  the  wing  whose  strength  corresponds  to  the  spanwise  variation 
of  the  circulation.  Using  a  linearized  model  that  neglects  roll-up,  we  assume 
that  the  vortex  sheet  coincides  with  the  wing  mean  surface,  defined  by 
Equation  (2c),  between  the  trailing  edge  and  downstream  infinity.  There  is  a 
jun^D  potential  function  across  the  sheet  which  is  constant  along 

lines  lying  in  the  sheet  that  are  parallel  to  the  freestream  direction.  The 
normal  component  of  velocity  and  the  pressure  must  be  continuous  through  the 
sheet,  however.  Thus,  on  the  vortex  sheet. 


)  (x  ,y  ,z  ^)  - 
te  o  m 


te  o  m 


r  (y  ),  4)  continuous, 
o  z 


(18) 


II 


Far  from  the  wing/body,  the  flow  is  undisturbed  except  in  the  downstream 
Trefftz  plane  where  the  vortex  sheet  induces  a  two-dimensional  downwash 
flow.  Thus, 

()>(“)  =  0  (19) 


everywhere  except  on  the  y-z  plane  at  downstream  infinity  where  the  potential 
function  satisfies  the  equation 


(a  ^  -  v^)  <ti  -  2vwi)  +(a^-w^)(J  =0 

00  '  ®  zz 


(20) 


subject  to  the  boundary  conditions 


(ji  -  (sin  o  +  )  = 

y  y  2 

surface 


=  0  . 


(21a) 


~  <ti(“»y»z”)  =  r(y)  on  the  Kutta  slit. 


(21b) 


^  >  0 


at 


/y^  + 


z2. 


(21c) 


The  Kutta  slit  corresponds  to  the  crossplane  profile  of  the  vortex  sheet  at 
infinity.  From  Equation  (15),  it  follows  that  a^  =  Boundary  condition 

(21a)  applies  on  a  cross-section  of  the  semi-inf initely  long  fuselage. 

For  the  unyawed  case,  only  half  of  the  configuration  need  be  considered, 
and  the  symmetry  condition  v  •  0  is  Imposed  on  the  vertical  plane  y  ■  0,  which 
contains  the  fuselage  centerline. 

Equations  (13)-(21)  collectively  define  the  problem  for  the  wlng/body 
flowfleld.  The  two-dimensional  problem  for  the  downwash  field  in  the  Trefftz 
plane  is  embedded  within  the  overall  three-dimensional  problem. 
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3 . 2  Computational  Domain 

The  computational  problem  is  formulated  by  transforming  the  equations  of 
Section  3.1  into  coordinates.  Substituting  Equations  (1)  and  (9)  into 

Equations  (13)  -  (14)  and  using  the  symbolic  definitions  given  in  Equation 
(10)  yields  the  equation  for  the  perturbation  potential: 


where 


*  (a^  -  G^C.a)  -  2„v  f«)  . 

c  (n) 


B  =  (a^  -  v^)  g(n)/b^  , 


(  2  2 
-  _  I  .  2  2.  7t2,j.  .  ,2  2,  'x2.-  .  ,  (a  -  w  ) 

I  T  (ri)c  (n) 


2uvG(C,n)G(E,n,i;)  -  |^uG(E.n)  +vG{^,n,4)j  h(0, 


D  =  I  j^(a^  -  v^)  G(C,n)  -  f(ng(n), 


o  r  9  7  —  VW  T 

=  I  (a  -  v^)  G(E,n,C)  -  uvG(E,n)  +  V(n)c(V)J 


F  =  2  (a^  -  u^)  +  (a^  -  v^)  G(C  ,n)G(C  ,n,C)  -  uv 


r((,n)G«.a)j  j^vGlt.a)  f«)h(C): 


13 
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Figure  5.  Schematic  diagram  of  computational  domain 
showing  boundary  conditions. 


L  = 


•Vy(C,ri)8(n)/b, 


G(C.n)  // (C,n)  +  G(?.n,C)  (C.r^) 


h(C), 


N  =  (cos  a)  -  sin  a. 


The  far-field  boundary  condition  remains  <t>  =  D  on  all  edges  of  the 
computational  domain  corresponding  to  infinity  except  the  transformed  Trefftz 
plane  where  the  governing  equation  becomes 


+  cT(h(l.  +  £^4*  =  JT, 

n  n  ^  ^  nc 


(25) 


with 


-  (af  -  v^)  g(n)/b^. 
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CT  = 


2  2-2  ~ 

(a„  -  V  )  G^l+C^,n,C)  ^ 


-  2vw 


T  (n)c  (n) 


G(l+C^,n,0 

T{n)c(n) 


h(0. 


(af  -  v^)  G(l-K  ,n,c)  -  .  .  -, 

“  o’  T(Ti)c(n) 


g(n)h(c). 


=  I 2vwl (n )  - 


(a^  -  v^)  H(14£ 


^.n,C)|  (w  -  sin  a)  T(n)c(n), 


and  velocities  v  and  w  follow  from  Equations  (23b)  and  (23c) 
with  (f)^  =  0.  Equation  (25)  is  subject  to  the  boundary  conditions 


(1-K  ,ri) 
o 


+ 


G(l-K^,n,c)  .^^(l-K^.n) 


1 

T(Ti)c(n) 


h(C)(|>^ 


-  sin  a  >  “  0, 

j  surface 


(26a) 


<ti  (C^g  ,n  ,0^)  -  h,0  )  =  r(Ti)  on  Kutta  slit, 


(26b) 


*  «  0  at  n  =  1+n  or  C  =  ±1*  (26c) 

o 

On  the  wing/body  symmetry  plane,  the  condition  v  =  0  is  enforced  using 
Equation  (23b). 


3.3  Local  Streamline  Coordinates 

For  use  in  applying  upwind-biased  finite  differences  (Section  4.1),  the 
potential  equation  is  rewritten  in  coordinates  that  are  locally  aligned  with 
the  stream  direction,  which  is  denoted  by  S.  In  such  coordinates,  the 
principal  part  of  Equation  (13)  takes  the  form; 

(a^  -  q^)  <t.gg  +  a^  (&4.  -  (tgg)  -  0,  (27) 
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where 
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T  (n  )c  (n ) 


V  G(C,n)  +  ' 


c<n) 


f(C)h{c), 


W  = 


2uvl (n)  +  V  H(r ,n ) 


u~  cosa  -  G(C,n)  (w  -  sin  a)  T(n)c(n) 


c(ri) 


+  <2v 


uI(C  ,n)  +  wl(n) 


+  H(C,ti,C);  (w  -  sin  a)  T(n)c(n). 


The  Laplacian  becomes 


where 


^  +G^(C,n) 


c^(n) 


f(?). 


“  g(n)/b  , 


"l  = 


G^(5,n)  +  G^a.n.c)  + 


T^(n)c^(n) 


hie), 


Sj  =  2G(C,n)f(f )8(n)/b. 


=  2G(^,n,i;)g(n)h(c)/b, 


Vj  =  2 


G(C.h) 

c(ri) 


+  G(F,n)  G(C,n,c) 


f(nh(c). 


(31) 
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Wj  =  H(f  ,n) 


cos  a  -  G(^,n)  (w  -  sin  a)  T(n)c(n)  c(n) 


+  (w  -  sin  a)  T(n)c(n). 


3 . 4  Transformation  Derivatives 

In  Sections  3.2  and  3.3,  various  quantities  appear  which  represent 
transformation  derivatives  between  the  (x,y,z)  and  (X,Y,Z)  coordinates.  These 
derivatives  follow  from  Equation  (1)  according  to  the  following  definitions: 


G(x,y ) 


H(x,y) 


Ky) 


G(x,y,z) 


11 

3y 


3^X 

3x3y 


3Z 


3y 


H(x,y,z) 


I(x,y) 


G(x,y ) 


3^Z 

3x3y 


11 

3x 


Ky) 


2 

3  Z 


3  y3  z 


G(C  ,n) 


H(£  ,n) 


I(n) 


G(C,n,c)  =G^(n,c)  +G2(f,n) 
H(f:  ,n,c)  =  Hj(n,c)  +  H2(Kb) 


I(C,n) 


G(C  ,n ) 


I(n) 


(32) 


(33) 


(34) 


(35) 


(36) 


(37) 


(38) 


(39) 
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Evaluation  of  these  transformation  derivatives  in  the  computational  domain 
uses  the  one-to-one  correspondence  that  exists  between  grid  points  in  the 
(x,y,z),  (X,Y,Z),  and  (f.  ,n,^)  coordinate  systems.  The  splitting  of  the 
functional  dependence  shown  in  Equations  <35)  and  (36)  actually  occurs  in  the 
(X,y,Z)  system.  Transformation  derivatives  which  are  not  Included  among 
Equations  (32)  -  (39)  either  are  zero  or  have  been  explicitly  evaluated  in  the 
equations  of  the  previous  sections. 
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4.0  NUMERICAL  SOLUTION  SCHEME 


4 . 1  Finite-Difference  Approximations 

The  type-dependent  finite-difference  concept  introduced  by  Murman  and 
Cole^^  is  the  basis  for  numerical  schemes  to  compute  steady-state  transonic 
flowfields.  Central  differences  are  used  to  approximate  the  potential 
equation  at  subsonic  points  of  the  solution  domain,  and  upwind-biased 
differences  are  used  at  supersonic  points.  Thus,  the  mathematical  character 
of  the  equation  is  properly  represented  as  it  changes  type  from  elliptic  to 
hyoerbolic.  The  original  application  involved  the  transonic  small-disturbance 
potential  equation. 

In  order  to  apply  upwind  differences  to  the  full  potential  equation,  it 
is  necessary  to  take  into  account  the  misalignment  between  coordinate  lines 
and  the  velocity  vector  at  any  given  point  of  the  flowfield.  The  principal 
part  of  the  equation  is  recast  in  a  form  that  constitutes  an  effective 
rotation  to  the  local  stream  direction,  denoted  by 

(I  -  +  (Ad>  -  (fgg)  =  0,  (40) 


where  M  =  q/a  is  the  local  Mach  number  and  the  quantities  ij)„  and  A(}>  are 
defined  in  Section  3.3.  Then,  upwind  differences  are  used  to  approximate 
contributions  to  the  first  term  on  the  lefthand  side  of  Equation  (40),  and 
central  differences  are  applied  to  factors  associated  with  the  second  term. 


Jameson  has  shown  further  that  the  relaxation  procedure  for  the  finite- 

difference  counterpart  of  Equation  (40)  can  be  viewed  in  terms  of  a  damped, 

three-dimensional  wave  equation  involving  an  artificial  time.^  The  time 

dependence  arises  because  of  the  appearance  in  difference  formulas  at  each 

grid  point  of  both  a  new  solution  value  from  the  current  relaxation 

/_  )  i  >  J  »  ^ 

step,  and  an  old  value  .  ,  from  the  previous  relaxation  step.  Thus, 

timelike  terms  occur  implicitly  which  correspond  to  These  terms  play  a 

role  in  controlling  the  stability  of  the  relaxation  process  but  do  not  affect 
the  final  solution  since  they  vanish  as  the  process  converges.  Sometimes  the 
damping  inherent  in  the  finite-difference  analog  of  Equation  (40)  must  be 


augmented;  this  can  be  accomplished  by  adding  to  Equation  (40)  a  term  of  the 
form 


‘ 


udt/u  V  .w  \ 

-  I  (S  +  —  (j>  +  O  1 

q  q  xt  q  yt  q  ztl 


(41) 


with 


f(C)u 

2  2 

q  c(n)A4 


P.  = 


c(n) 


+  vG(4  ,n) 


,  -  f(C)u  V 

^2  “  2  ,  b  ’ 

q  c(ri)AC 


"2  =- 


ilSju 


q  c(ri)A4 


uG(n)  +  vG(C,n,4)  + 


w 


T  (n)c(n) 


In  Equation  (41),  the  value  of  parameter  £  Is  arbitrary  and  can  be  adjusted  to 
control  the  amount  of  damping  augmentation. 


In  writing  difference  approximations,  central  differences  based  on  old 
values  are  used  to  calculate  first  derivatives  required  to  evaluate  the 
velocity  components  in  Equations  (23a)  -  (23c).  At  grid  points  where  the  flow 
is  subsonic.  Equation  (22)  is  used,  and  second  derivatives  are  represented  by 
central  differences.  Typical  formulas  are 
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(fv,),  = 


_ .f  \.(n) 

2(^(r  )2  i+l,j,k  i,J,k  j 


-  *  "'i.j.k  [s  ‘["Ilk  ^  (‘  -  =) 

(n+l)  ) 


i.j.k  i-l,j,k  I 


(42) 


1 


2(Ari)^ 


(g  +  g  )  (rf)  ) 


(g  +  g  )  ) 


(43) 


U 


1 

4ACA<; 


(n) 

i+1 , j ,k+l 


(n)  ^^(n+1) 

''^i+l,j,k-l  ^i-l,j,k-l 


.(n+l) 

'^i-1 ,  j  ,k+l 


(44) 


In  Equation  (42),  the  relaxation  factor  w  has  been  incorporated  into  the 
difference  expression.  Also,  in  Equations  (42)  and  (43),  stretching  function 
values  required  at  half-step  points  have  been  evaluated  as  the  average  of  the 
values  at  grid  points  on  either  side;  the  averaging  is  a  second-order 
approximation  to  the  stretching  function  at  the  half-step  location.  At  grid 
points  where  evaluation  of  the  local  velocity  indicates  the  flow  to  be 
supersonic.  Equation  (40)  is  used,  and  second  derivatives  contributing  to 
H)gg  in  the  first  term  are  upwind  differenced.  For  example,  in  the  case  where 
velocity  components  u,  v,  and  w  are  all  positive. 


(f.,),  = 


2(AO 


^U,j,k  "i-l,j,k^  '  ^i.j.k  ^l,j,k  ^1-1, j.k^ 


1-1, j,k  i-2,j,k  i-l.j,k  i-2,j,k  J 


(45) 
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A4AC  ['^i.J.k  ~  '*^i,j,k-l  '*'  ^i-l,j,k-l  “  '^l-l,j,k 


(46) 


Derivatives  associated  with  the  ijig^  terms  are  approximated  by  upwind 
differences;  a  representative  form  is 


1-1. j.k  L 


_  (^(^+1)  -  Z*')  )1 


(47) 


If  u,  V,  or  w  is  negative,  Equations  (45)  -  (47)  are  revised  to  reflect  the 
appropriate  upwind  direction.  All  second-derivative  contributions  to  the 
term  (Aif  -  4>  )  in  Equation  (41)  are  central  differenced  in  a  manner  similar 

Od 

to  Equations  (43)  and  (44). 

4.2  Boundary  Conditions 

Since  the  present  formulation  produces  a  finite  computational  domain 
(Figure  5),  application  of  far-field  conditions  is  straight-forward  on  those 
faces  that  represent  infinity  and  on  which  <p  ••  0»  The  boundary  condition  in 
the  Trefftz  plane  (downstream  infinity  face)  depends  on  the  circulation 
distribution  and  is  not  known  a  priori;  therefore,  it  must  be  calculated  by  an 
embedded  relaxation  procedure  which  solves  the  two-dimensional  problem  for  the 
downwash  field.  This  calculation  is  based  on  Equation  (25)  with  boundary 
conditions  (26a)  -  (26c) ;  central-difference  approximations  are  used  similar 
to  those  at  subsonic  grid  points.  The  Tref f tz-plane  calculation  coincides 
with  updates  of  the  circulation  distribution.  Finite  differences  in  the  4- 
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direction  that  traverse  the  vortex  sheet  are  adjusted  to  incorporate  the 
potential  junqj  across  it.  At  the  wlng/body  symmetry  plane,  where  v  =  0, 
Equation  (23b)  is  used  to  compute  values  of  the  potential  function  at  image 
points  located  one  grid-step  beyond  the  computational  domain.  In  order  to 
compensate  for  coordinate  nonorthogonality  at  this  plane,  derivatives  are 
replaced  by  one-sided  differences  that  are  oriented  to  correspond  with  the 
skewed  intersection  of  spanwlse  grid  lines.  The  present  symmetry-plane  scheme 
is  an  extension  of  one  used  by  Boppe^  and  has  proved  to  be  numerically  stable. 

At  the  wlng/body  surface,  special  treatment  is  required  to  enforce  the 
flow  tangency  condition  given  by  Equation  (24).  A  scheme  is  used  that  is 
based  on  ideas  derived  from  previous  numerical  applications  of  non-surface- 
fitted  coordinates . The  Neumann  form  in  Equation  (24)  is  replaced  by  an 
equivalent  Dlrichlet  condition  that  is  applied  at  image  points  below  the 
wing/body  surface.  Referring  to  Figure  6,  point  S  is  a  typical  surface 
control  point  defined  by  the  intersection  of  a  vertical  grid  line.  Using  the 
corner-shaped  boundary-point  array  shown,  one-sided  differences  based  on  the 
Lagrange  interpolation  formula  for  three  unevenly  spaced  points  are 
substituted  into  Equation  (24)  to  obtain  a  Dlrichlet  value  for  the  potential 
function  at  point  S.  Extrapolation  along  line  S-1-2  then  transfers  this  value 
to  the  uniformly  spaced  image  point  D1  located  below  the  surface  boundary. 

Thus , 

<(i(Dl)  =.>'pA(S),  -r(S),  A5,  An,  AC,  aJ  ,  (48) 

where //^(S),  //^(S)  are  known  surface  slopes,  are  potential  values  at 

indicated  points  of  the  boundary  array  associated  with  point  S,  (AE,An,AC)  are 
grid  stepsizes  in  the  respective  coordinate  directions,  and  6  is  the  offset 
distance  between  point  S  and  the  grid  point  above  it.  Repeated  applications 
of  the  Lagrange  formula  to  the  sets  of  points  connected  by  arrows  in  Figure  6 
provide  the  inter-grid  values  4i  . .  .4)  .  The  potential-function  value  at  Image 
point  D1  fixes  the  boundary  condition  during  relaxation  of  the  solution  at 
exterior  points  on  the  vertical  line  DI-2.  After  each  relaxation  sweep  of  the 
exterior  field,  the  image-point  potential  value  is  recomputed.  In  the  event 
that  upwind  differencing  at  a  surface-adjacent  point  necessitates  a  second 
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5  Surface  point  at  which  Neuman  flow-tangency 
condition  is  enforced 

1*6  Exterior  points  used  to  calculate  equivalent 
Oirichlet  potential  at  surface  point  S 

D1,  02  Interior  image  points  used  in  numerical 
solution  scheme 

6  Grid  point  offset  distance 


QPO3>110O-e 


Figure  6.  Grid>point  array  used  to  numerically  enforce  the 
boundary  condition  at  a  typical  wing/body  surface  point. 


image  point  D2,  the  associated  potential  value  is  obtained  by  linear 
extrapolation  in  the  vertical  direction. 

In  order  to  treat  each  control  point  consistently,  the  corresponding 
boundary-point  array  is  usuf ily  directed  to  the  exterior  side  of  the  wing/body 
surface.  Thus,  a  reorientation  of  the  array  occurs  at  a  locus  of  zero 
streamwise  slope  on  a  convex  part  of  the  surface,  as  shown  in  Figure  7.  This 
arrangement  has  the  effect  of  eliminating  direct  coupling  between  neighboring 
image  points  in  the  present  methodology.  The  boundary-point  array  is  not 
reoriented  in  concave  surface  regions,  however. 

Additionally,  it  has  been  found  useful  to  permit  surface  image  points  to 
be  multi-valued  in  order  to  more  accurately  represent  the  surface  condition  at 
configuration  edges.  Figure  8  indicates  schematically  the  situation  in  a 
wing-section  plane.  Image  point  P  represents  the  leading-edge  condition 
during  relaxation  of  the  solution  on  the  line  ILE-1  upstream  of  the  section. 
Then,  for  computation  on  the  segment  of  line  ILE  below  the  section,  point  P 
represents  a  lower-surface  condition.  In  the  same  way,  point  P  can  also  be 
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assigned  an  upper-surface  boundary  value.  Image  point  Q  at  the  section 
trailing  edge  is  treated  similarly.  Other  such  multi-valued  image  points 
occur  at  the  fuselage  nose,  along  the  side-edge  locus  of  the  fuselage,  and  at 
the  wing  tip. 

4.3  Computation  Procedure 

The  system  of  algebraic  difference  equations  for  the  wing/body  problem  is 
solved  iteratively  by  systemlcally  sweeping  through  the  computational  space 
(Figure  5)  and  relaxing  the  solution  along  vertical  columns  of  grid  points. 

In  the  present  methodology,  the  domain  is  swept  by  crossplanes  beginning  at 
the  upstream-infinity  boundary  and  proceeding  to  the  downstream-infinity 
boundary.  Each  crossplane  is  swept  by  vertical  lines  from  the  configuration 
symmetry  plane  to  the  spanwise-lnf Inity  face  of  the  domain.  In  a  crossplane 
that  Intersects  the  wing/body  surface,  all  column  segments  beneath  the 
configuration  are  relaxed  first,  followed  by  column  segments  above  the 
geometry,  and  then  by  columns  in  the  outboard  region  of  the  crossplane.  The 
relaxation  procedure  is  continued  until  either  a  prescribed  convergence 
criterion  is  satisfied  or  a  specified  number  of  domain  sweeps  have  been 
completed. 

During  a  particular  iterative  cycle,  the  first  step  is  to  fix  the 
surface-boundary  condition  by  calculating  image-point  potential  values 
associated  with  each  surface  control  point.  The  solution  is  then  relaxed 
throughout  the  domain  as  described  above.  At  designated  intervals,  the 
circulation  distribution  is  updated  by  applying  Equation  (26b)  at  the  wing 
trailing  edge.  Following  each  circulation  update,  the  Tref f tz-plane  boundary 
condition  is  recalculated  by  an  embedded  relaxation  of  Equations  (25)  -  (26) 
for  the  downwash  field.  The  cycle  is  then  repeated. 

During  the  sweep  process,  image-point  values  for  the  wlng/body  surface 
condition  are  substituted  sequentially,  as  required,  into  the  solution  array 
that  stores  the  potential  function.  This  procedure  simplifies  program  logic 
by  effectively  eliminating  the  distinction  between  surface-adjacent  grid 
points  and  interior  field  points  in  computing  finite  differences.  On  vertical 
grid  lines  that  Intersect  the  wing/body  configuration,  grid  benchmarks  for  the 
surface-adjacent  points  are  used  to  designate  the  length  of  column-segments  on 
which  the  solution  is  relaxed. 
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5.0  RESULTS 


Representative  calculations  have  been  made  for  a  configuration  comprised 
of  ONERA  Wing  M6  attached  at  mid-height  to  a  hemisphere-cylinder  fuselage. 
Details  of  the  wing  geometry  are  given  in  References  18  and  19.  The  wing  has 
a  planform  with  a  leading-edge  sweep  angle  of  30°,  a  taper  ratio  of  0.56,  and 
a  uniform  section  whose  thickness  ratio  is  0.098.  The  fuselage  radius  is 
taken  to  be  25%  of  the  exposed  wing  semispan.  A  representation  of  the 
configuration  planform  is  shown  in  Figures  9  and  10. 

Figures  9  and  10  also  show  the  calculated  pressure  distributions  at 
several  spanwise  stations  on  the  configuration  for  a  freestream  Mach  number  of 
0.84  and  angles  of  attack  of  0°  and  3.06°,  respectively.  The  sharp  peaks  in 
the  wing  leading-edge  region  are  the  consequence  of  grid  coarseness  (see 
below).  There  is  an  indication  of  a  double-shock  structure  in  the  midspan 


Figure  9.  Calculated  surface-pressure  distributions  on  a  configuration  composed  of  ONERA  Wing  M6 
attached  at  mid-height  to  a  hemisphere-cylinder  fuselage:  =  0.84,  a  =  0°.  Grid  dimension: 

49  x  19  X  25.  Every  third  spanwise  station  is  shown. 
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Figure  10.  Calculated  surface-pressure  distributions  on  a  configuration  composed  of  ONERA  >^ing 
M6  attached  at  mid-height  to  a  hemisphere-cylinder  fuselage:  =  0.84,  a  =  3.06°.  Grid  dimension: 

49  X  19  X  25.  Every  third  spanwise  station  is  shown. 

region  of  the  wing,  which  coalesces  to  a  single  strong  shock  near  the  wing 
tip;  this  effect  is  particularly  evident  for  the  lifting  case  in  Figure  10. 

In  both  cases,  the  pressure  distribution  at  the  fuselage  centerline  is 
modified  by  the  presence  of  the  wing,  again  with  a  more  pronounced  effect  in 
the  lifting  case. 

In  order  to  assess  solution  accuracy,  the  calculations  are  compared  with 
available  experimental  data  for  ONERA  Wing  M6  obtained  in  a  wing-alone  test 
(References  18  and  19).  Since  spanwise  grid  stations  and  test  stations  on  the 
wing  do  not  coincide,  the  computed  results  are  interpolated  linearly  along 
constant-percentage-chord  lines  to  the  positions  of  the  data  measurements. 

The  comparisons  are  shown  in  Figures  11  and  12.  In  general,  agreement  is 
reasonably  good  especially  near  the  wingtlp  where  the  presence  of  the  fuselage 
in  the  calculations  has  the  least  effect.  Differences  near  the  wing  leading 
edge  and  the  poor  resolution  of  the  double-shock  structure  along  the  wing  can 
be  attributed  to  grid  coarseness  in  the  computed  results,  while  traillng-edge 
discrepancies  are  more  probably  the  consequence  of  viscous  effects  in  the 
data.  (Note  the  definite  shock-induced  separation  exhibited  by  the  test  data 
at  the  wingtlp  station  in  Figure  12.)  Anomalous  behavior  such  as  the 


30 


The  calculations  for  these  cases  were  carried  out  on  a  49  x  19  x  25  grid 
in  the  chordwise,  spanwise,  and  vertical  directions,  respectively.  Thirty™ 
three  vertical  grid  lines  intersect  each  wing  section  chord,  thus  providing 
sixty-six  (upper  plus  lower)  surface  control  points  on  each  section  profile. 
Control  points  are  spaced  at  approximately  every  3%  of  section  chord  with  an 
edge-offset  of  0.5Z  chord  at  the  leading-  and  trailing-edge  points.  The 
pressure  peaks  in  Figures  9  and  10  correspond  to  the  second  chordwise  control 
point  from  the  leading  edge  on  each  section.  In  the  spanwise  direction,  13 
grid  stations  occur  on  the  configuration  semispan. 

The  calculations  were  performed  on  a  Control  Data  CYBER  175  computer  with 
FTN(0PT=2)  compiler.  The  nonlifting  case  of  Figure  9  required  2.10  min  of  CPU 
execution  time,  and  the  lifting  case  of  Figure  10  required  5.23  min. 
Convergence  is  based  on  a  10”^  cut-off  limit  on  the  maximum  correction  to  the 
potential  function  over  the  solution  domain  between  the  final  two  iterative 
sweeps.  Decreasing  the  cut-off  limit  to  10“^  was  found  to  approximately 
double  solution  times. 

Table  1  shows  the  accuracy  of  the  scheme  used  to  numerically  apply  the 
surface  boundary  condition.  At  two  semispan  positions  on  the  wing,  velocity 
slopes  computed  from  a  converged  solution  are  compared  with  prescribed 
streamwise  surface  slopes.  Agreement  between  corresponding  slope  values 
extends  to  at  least  the  third  decimal  place,  one  order-of-magni tude  greater 
than  the  cut-off  limit  on  the  potential  function. 

In  order  to  explore  the  stability  of  the  computer  program,  additional 
calculations  were  performed  for  a  configuration  with  the  planform  shown  in 
Figure  A  in  which  the  wing  has  a  uniform  NACA  0012  section.  Convergent 
operation  without  the  use  of  damping  augmentation  was  achieved  in  nonlifting 
cases  for  freestream  Mach  numbers  up  to  0.99  and  at  selected  (supercritical) 
Mach  numbers  for  angles  of  attack  up  to  6®.  These  results  are  not  included  in 
this  report  since  no  comparison  data  exist  for  the  particular  geometry. 


*■  of  surface-slope  values  at  two  semispan  stations. 
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2Y/8  - 

.2457 

2Y/B  • 

.9900 

tK-XLEHC 

SLOPE  U 

OZOKU 

SLOPE  U 

OZDXU 

.0050 

1.191132 

1.192676 

1.192739 

1.192676 

*  0  3  6  S 

. 199445 

.199572 

.199503 

.199572 

.  06  7 7 

. 137938 

.133035 

.136046 

.138035 

.0^89 

.109128 

.109208 

.109222 

.109208 

•  L  3  0  0 

.087512 

.087580 

.08  7594 

.037580 

•  1611 

,071034 

.071 092 

.071107 

.071092 

.1971 

.05  8087 

.058136 

.058151 

.050136 

.  2730 

.047324 

,047365 

.047380 

.047365 

.  2539 

.037729 

.037764 

.03  7777 

.037764 

.  23S8 

.028551 

.028580 

.028593 

.028530 

•  3 1  5  6 

.019317 

.019339 

.019351 

.019339 

.  3  ’» 6  4 

.009725 

.009742 

.009751 

.009742 

.  3771 

-.000267 

-.000277 

-.000303 

-.000277 

.4079 

-.010686 

-.010661 

-.010620 

-.010661 

.4  386 

-.02  1286 

-.021235 

-.021118 

-.021235 

.  4693 

-.031790 

-.031724 

-.03  1526 

-.031724 

.5000 

-.041930 

-.041359 

-.04  15  78 

-.041859 

.  5  3  0  7 

-.051393 

-.051328 

-.050965 

-.051329 

.5614 

-.059962 

-.0599;  :i 

-.0594  73 

-.059911 

.  5971 

-.0674  39 

-.067410 

-.066905 

-.067410 

.6229 

-.073775 

-.073773 

-.07  32  10 

-.073773 

.  65  36 

-.079030 

-.079056 

-.078448 

-.079056 

•  6  6  4  4 

-.033393 

-.063447 

-.032803 

-.083447 

.7152 

-.037183 

-.097268 

-.036599 

-.037269 

.  7461 

-.090812 

-.090916 

-.090227 

-.090713 

.  7770 

-.094731 

-.094861 

-.094151 

-.094861 

.  8079 

-.099350 

-.099501 

-.098771 

-.099501 

.  8309 

-.104921 

-.106092 

-.104338 

-.105092 

.  8700 

-.111292 

-.111480 

-.110698 

-.111490 

.  9011 

-.117739 

-.117940 

-.117130 

-.117940 

.9  32  3 

-.122597 

-.122904 

-.121977 

-.122804 

.  963  6 

-.123504 

-.123781 

-.122973 

-.123781 

.9950 

-.12  3582 

-.123776 

-.122919 

-.123776 

(dij  Convercience  lirrur  on  ^0;  iQ  ^ 

(b)  Nonitfting  C3lcul3tion  for  ONERA  Wing  M6/h6nTisph©re-cvlio(jpr  fus^idge 


34 


6.0  CONCLUDING  REMARKS 


The  present  work  substantiates  the  use  of  non-surf ace-f itted  coordinates 
for  numerical  computation  of  transonic  wing/body  flowfields.  However,  a 
number  of  modifications  and  improvements  remain  to  be  incorporated  in  order  to 
make  the  present  computer  program  useful  for  engineering  applications. 

The  capability  of  the  program  to  perform  calculations  on  a  series  of 
progressively  finer  grids  needs  to  be  completed  in  order  to  improve  computing 
efficiency.  This  modification  involves  two  requirements.  Verification  of  the 
grid-halving  subroutine  in  the  program  must  be  completed,  and  the  array  that 
stores  the  potential  function  must  be  moved  to  disk  storage,  with  a  sequential 
transfer  into  central  memory  of  only  those  array  segments  required  at  each 
stage  of  the  computation  process. 

The  chordwise  stretching  function  used  between  the  leading  and  trailing 
edges  of  the  wing  is  symmetric  about  the  mid-chord  locus.  In  order  to  Improve 
resolution  of  the  blunt  leading-edge  region  within  a  fixed  grid  dimension,  the 
Introduction  of  an  asymmetric  stretching  that  clusters  more  points  near  the 
leading  edge  than  the  trailing  edge  would  be  helpful.  This  modification  would 
provide  a  more  efficient  means  of  improving  solution  accuracy  in  the  leading- 
edge  region  than  a  simple  Increase  in  the  number  of  chordwise  grid  points. 

Incorporation  into  the  program  of  a  finite-fuselage  capability  would 
provide  a  better  geometry  model  for  engineering  applications.  An  attempt  to  do 
this  during  the  present  contract  was  unsuccessful  when  the  resulting  version 
of  the  code  proved  to  be  nonconvergent .  The  principal  difficulty  seemed  to 
involve  the  question  of  how  to  properly  treat  the  vortex  sheet  (circulation 
distribution)  in  the  region  behind  the  fuselage. 

Finally,  there  exists  the  problem  of  improving  flowfleld  and  geometry 
resolution  around  the  fuselage.  A  method  for  accomplishing  this,  depicted  in 
Figure  13,  involves  a  two-coordinate  formulation  in  which  the  present 
coordinate  arrangement  is  retained  about  the  wing  and  better-suited 
coordinates  —  perhaps,  a  simple  stretched,  cartesian  system  —  are  introduced 
in  a  vertical  slab  whose  width  coincides  with  that  of  the  fuselage.  This 
scheme  would  require  interpolation  within  an  overlap  region  (possibly,  only  a 
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single  plane)  common  to  both  coordinate  systems-  In  order  to  conveniently 
accommodate  the  two-coordlnate  formulation.  It  may  be  preferable  to  revise  the 
computation  strategy  such  that  the  domain  is  swept  by  wing-section  planes 
beginning  at  spanwise  infinity  and  moving  to  the  wing/body  centerline.  The 
vertical-line  sweep  in  each  section  plane  would  begin  at  the  upstream  boundary 
and  proceed  to  the  downstream  boundary. 

The  baseline  computer  program  described  in  Appendix  A,  when  modified  as 
discussed  above,  is  expected  to  provide  a  framework  for  treatment  of  complex 
wing/body  configurations.  Further  extension  to  include  add-on  components  such 
as  nacelles,  stores,  or  additional  lifting  surfaces  should  be  possible  also. 


GP  »  T  05^7  8 

Figure  13.  Two-grid  arrangemenl  for  wing/body  flow 
analysis. 
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APPENDIX  A.  USER'S  GUIDE  TO  COMPUTER  PROGRAM 


A  computer  program  based  on  the  formulation  presented  in  the  main  body  of 
this  report  is  listed  in  Appendix  B.  The  subroutine  structure  of  the  program 
is  shown  in  Figure  Al,  and  a  summary  of  subroutine  functions  is  given  in  Table 
Al.  The  code  has  been  run  on  a  Control  Data  CYBER  175  computer  with  an  FTN 
compiler.  As  dimensioned,  it  requires  260k  (octal)  storage  locations  to  load 
and  execute.  The  present  form  of  the  program  does  not  use  peripheral  storage 
devices. 

The  solutions  presented  in  this  report  have  been  computed  on  a  single 
grid.  Some  logic  is  contained  in  the  code  to  permit  a  calculation  to  be  made 
on  a  series  of  progressively  finer  grids,  but  it  has  not  been  fully 
Implemented.  In  particular,  subroutine  HALFS  interpolates  the  potential  array 
P(1,J,K)  from  an  initial  grid  onto  one  that  is  half-spaced  in  each  coordinate 
direction  of  the  computational  domain.  The  interpolated  array  is  then  used  as 
the  starting  solution  for  continuation  of  the  relaxation  process  on  the  new 
grid.  Input  parameter  NHALF  specifies  the  number  of  grid-halving  cycles  to  be 
performed.  Full  implementation  of  this  capability  will  require  that  the  array 
P(I,J,K)  be  transferred  to  disk  storage  and  that  provision  be  made  for  a 
sequential,  plane-by-plane  transfer  of  P(I,J,K)  values  into  central  memory  to 
coincide  with  each  relaxation  sweep  through  the  computational  domain. 

Figure  1  defines  the  class  of  wlng/body  geometries  that  can  be 
represented  by  the  program.  Input  data  are  smoothed  to  ensure  that  the  wing 
leading  and  trailing  edges  are  piecewise  straight  lines.  One  break  (kink)  is 
permitted  in  each  edge  but  need  not  be  present.  If  both  leading-  and 
trailing-edge  breaks  occur,  they  may  be  at  different  spanwlse  positions. 
Although  the  formulation  outlined  in  Section  2.0  places  no  restriction  on  wing 
attachment  to  the  fuselage,  code  logic  assumes  that  the  horizontal  mid-plane 
of  the  grid  (5  -  0)  coincides  with  the  maximum  width  position  on  the 
fuselage.  Combined  with  the  specification  of  circular  fuselage  cross- 
sections,  this  computational  arrangement  effectively  limits  application  to 
mid-wing  configurations  at  present.  This  restriction  can  be  removed  by 
incorporating  a  two-coordinate  formulation  as  discussed  in  Section  5.0. 
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Figure  AI.  Subroutine  structure  of  the  wing/body  program. 
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TABLE  Al.  SUBPROGRAM  LIST 


NAME 

FUNCTION 

MAIN 

Reads  input  data;  controis  overaH  program  logic 

HALFS 

Interpolates  solution  arrays  onto  a  half-spaced  grid.  (Has  not  been  verified.) 

GRID 

Defines  computational  grid;  calculates  wing  configuration  data  at  spanwise  grid  stations 

GEOM 

Calculates  wing/body  coordinates  and  slopes  at  planform  grid  stations  in  the  physical  domain 

ASPLINE 

Parameterizes  input  geometry  data  in  terms  of  arc  length  along  the  curve  prior  to  spline  interpolation 

SSPLINE 

Interpolates  input  geometry  data 

PROFL 

Calculates  geometry-defining  quantities  in  the  computational  domain 

COEFF 

Calculates  fixed  quantities  which  appear  in  the  finite-difference  equations 

(Nir 

Initializes  solution  arrays 

SOLVE 

Executes  one  relaxation  sweep  of  the  computation  domair  ;  updates  the  circulation  distribution  and 
solution  boundary  values 

FARBC 

Computes  the  Treffiz-plane  boundary  condition 

SURFBC 

Calculates  boundary  values  at  control  points  on  the  wmg/body  surface 

SURFSC 

(Entry)  calculates  boundary  values  at  side-edge  points  of  the  configuration 

SYMBC 

(Entry)  Calculates  image  values  at  the  wmg/body  symmetry  plane 

TRICOE 

Calculates  coefficients  of  the  f mite-difference  potential  equation  at  grid  points  along  a  specified  vertical 
grid-line  segment 

TRIT 

(Entry)  calculates  coefficients  of  the  finite-difference  downwash  equation  m  the  Trefftz  plane 

INVERT 

Solves  the  finite-difference  potential  equation  along  a  vertical  gnd-lme  segment  to  obtain  updated  solution  values 

MAXI 

Determines  the  maximum  potential-value  increment  along  a  specified  vertical  gnd-line  segment  between  the 

Current  and  preceding  relaxation  sweeps 

ARRAYS 

Prints  out  solution  arrays  to  specified  iteration  intervals:  circulation  distribution,  potential  distribution, 

Surface  boundary  values  and  associated  image-point  values,  symmetry-plane  image-pomt  values.  (Used  mainly 
for  diagnostic  purposes.) 

CPCOMP 

Calculates  and  writes  the  pressure-coefficient  distribution  on  the  wmg/body  surface 

QP0JH901 


Table  A2  summarizes  the  sequence  and  format  of  input  data  required  by  the 
program.  Data  categories  are  as  follows: 

•  case  title, 

•  computational  grid  parameters, 

•  code  execution  parameters, 

•  case  specification  (Mach  number,  angle  of  attack), 

•  fuselage  configuration  data,  and 

•  wing  configuration  data. 

Suggested  values  for  grid  and  execution  parameters  are  included  in  the 
table.  Also  given  are  parameter-value  limits  imposed  by  current  code 
dimensions . 
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TABLE  A2.  GLOSSARY  OF  INPUT  DATA 


CARO 

COLUMNS 

VARIABLE 

EXPLANATION 

1 

1  80 

TITLE 

Case  title  (written  to  output) 

2 

1-10 

NXI 

Number  of  chordwise  grid  steps  at  start  of  calculation.  Maximum:  48 

11-20 

A1 

X  —  ^  stretching  constant  in  £q.  (3b).  Suggested  value;  0.16 

21-30 

A2 

X  —  ^  stretching  constant  m  Eq.  (3bl.  Suggested  value:  2  75 

31-40 

XIO 

Stretching  transition  point  of  ^  ~  coordinate;  see  Figure  2.  Ratio  XIO/ (1  +  XIO) 
determines  fraction  of  chordwise  grid  steps  which  occur  on  each  wng-section  chord 

41  50 

XCAPO 

Stretching  transition  point  of  X-coordinate;  see  Figure  2.  Must  be  inside  wing-section 
edge.  Suggested  value.  0.495 

3 

1  10 

Number  of  spanwtse  grid  steps  at  start  of  calculation.  Maximum;  1 8 

11  20 

21  30 

Y  —  T?  stretching  constant  in  Eq.  (5b).  Suggested  value;  0.16 

Y  -  T? stretching  constant  in  Eq.  (5b).  Suggested  value:  2.75 

31  40 

ETAO 

Stretching  iiansition  point  of  T?  —  coordinate,  see  Figure  2.  Ratio  ETAOd  +  ETAO) 
determines  fraction  of  spanwise  grid  steps  which  occur  between  the  wing/body 
symmetry  plane  and  the  wmgtip 

41  60 

YCAPO 

Stretching  transition  point  of  Y  coordmate;  see  Figure  2.  Must  be  inside  wmgtip. 

Suggested  value .  0.49995 

4 

110 

N2ETA 

Number  of  grid  steps  m  vertical  direction  at  start  of  calculation.  Maximum:  24 

\t  20 

Cl 

Z  —  f  stretching  constant  in  Eq.  (7)  Suggested  value;  0.45 

5 

1  10 

(TERM 

Maximum  number  of  iterations  to  be  executed  or^  the  initial  grid 

11  20 

NHALF 

Number  of  grid  halving  cycles.  Set  NHALF  0  for  single  grid  calculation,  (Note; 

Subroutine  HALFS  has  not  been  fully  verified  1 

21  30 

NPRINT 

Iteration  frequency  for  execution  of  subroutine  ARRAYS,  which  prints  complete 
solution  arrays  for  diagnostic  purposes.  If  NPRINT  ■  0.  only  the  circulation  distribution 

IS  printed  open  cofnp/etior  of  the  feloMation  procedure  If  NPRINT  >ITERM,  complete 
solution  arrays  are  printed  after  the  relaxation  process 

6 

1  10 

WE 

Relaxatioi'  parameter  for  potentia*  a»  elliptic  field  pomts.  Suggested  value.  1.70 

11  20 

WG 

Relaxation  parameter  for  Circulation.  Suggested  value;  1.00 

21  30 

DPLIM 

Convergence  cut  off  limit  for  maximum  potential-value  increment  between  successive 
Iterations.  Suggested  value  lo*^ 

31  40 

EPSI 

Damping  factor  m  potential  difference  eriuation  used  at  hyperbolic  field  points. 

Suggested  value  0.00,  increase  if  instability  occurs 

7 

MO 

2MACH 

Freestream  Mach  number 

II  20 

ALPHA 

Angle  of  attack  ot  the  wmg  reference  plane  (m  degrees) 

s 

1  10 

NF 

Number  of  fuselage  coordinate  sets  to  be  read  from  the  following  cards  (one  XF,  ZF, 

RF  set  per  card) 

INF 

1  10 

XF 

X  coordinate  along  fuselage  beginning  at  nose,  see  Figure  1 

11  20 

ZF 

/  coordinate  of  fuseiage  side-profile  reference  line;  see  Figure  1 

21  30 

RF 

Fuselage  ctossplane  radius;  see  Figure  1 

9 

1  10 

JSECT 

Number  of  wing-seciion  data  sets  to  be  read  subsequently.  Maximum  5 

11  20 

JBL 

Sequence  number  of  wing  data  set  at  break  in  leading  edge.  If  no  break  exists,  set  JBL  -  0. 

21  30 

JBT 

Sequence  number  ot  wmg  data  set  at  break  m  trailing  edge.  It  no  break  exists,  set  JBT  0. 

The  fust  rdfci  uses  .ilphanumenc  format  20A4.  All  remaining  cards  use 
repeated  f loatmq  point  format  8F  10  5  Cor^voision  of  data  to  irueger 
mode  IS  perfomrecJ  within  the  program  as  required 
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TABLE  A2.  (CONTINUED)  GLOSSARY  OF  INPUT  DATA 


CARO 

COLUMNS 

VARIABLE 

10 

1-10 

YS 

XLES 

ZLES 

CS 

41-50 

ATS 

51-60 

TS 

61-70 

FS 

11 

1-10 

NSU 

11-20 

NSL 

21-30 

i 

KSYM 

i 

1-NSU 

1-10 

XSU 

11-20 

ZSU 

1-NSL 

1-10 

XSL 

11-20 

ZSL 

12 

- 

- 

13 

- 

- 

U-IS") 

16-17  y 

18-19  J 

- 

EXPLANATION 


Spanwise  station  ot  wing  section 
x-coordinate  of  wing-section  leading  edge 
2-coordinate  of  wing-section  leading  edge 
Wing-section  chord  length 

Wing-section  twist  angle  relative  to  wing  reference  plane  (in  degrees) 

Wing-section  thickness-to-chord  ratio 

Repeat  indicator.  If  FS  =  0,  wing-section  coordinate  data  from  the  previous  span  station 
IS  used:  the  next  card  specifies  parameters  for  the  wing  section  at  the  next  span  station 
(go  to  Card  12).  If  FS  =  1,  coordinates  for  a  new  section  profile  are  read  from  the 
following  data  cards. 

Number  of  wing-section  upper-surface  coordinate  sets  to  be  read  from  the  following 
cards  (one  XSU,  ZSU  set  per  card) 

Number  of  wing-section  lower-surface  coordinate  sets  to  be  read  from  the  following 
cards  (one  XSL,  2SL  set  per  card) 

Wing-section  symmetry  indicator.  If  KSYM  =  0,  the  section  is  asymmetric;  both  upper- 
and  lower-surface  coordinates  must  be  given.  If  KSYM  =  1 ,  the  section  is  symmetric,  and 
only  the  upper-surface  coordinates  are  required. 

Upper-surface  coordinates  of  wing  section  from  leading  to  trailing  edge.  One  set 
per  card 

Lower-surface  coordinates  of  wing  section  from  leading  to  trailing  edge.  One  set  per  card. 
(Required  only  if  KSYM  =  0) 

Repeat  of  Data  Card  10  for  the  next  wing  section 

Repeat  of  Data  Card  1 1  and  data  sets  f  -NSU  and  1  -NSL  for  the  wing  section  defined  by 
Data  Card  12.  (Required  only  if  FS  =  1  on  Card  1 2) 

Up  to  three  additional  wing-section  definition  cards  and  coordinate  data  sets.  Total  number 
of  data  sets  must  correspond  to  JSECT  on  Data  Card  9. 


Fuselage  data  consist  of  sets  of  side-profile  reference-line  coordinates 
and  crossplane  radii  given  at  strearawlse  stations  beginning  at  the  nose  and 
proceding  aft.  These  data  are  spline  interpolated  to  obtain  required  values 
at  grid-point  stations.  Along  a  constant-radius  segment  of  the  fuselage,  a 
number  of  values  must  be  given  to  ensure  accuracy  of  the  spline  fit,  and  the 
data  must  be  more  closely  spaced  near  the  ends  of  the  segment  than  in  its 
central  region.  A  final  table  value  at  streamwise  location  XF  >  25  is 
required  to  specify  the  semi-infinite  fuselage  length. 

Wing  data  are  read  sect lon-by-sect ion  beginning  at  the  wing  root  and 
moving  outboard  to  the  wingtlp.  In  the  data  set  for  each  section,  the  first 
card  specifies  section  properties:  spanwise  position,  leading-edge 
coordinates,  chord  length,  twist  angle  relative  to  a  horizontal  reference 


line,  thickness-to-chord  ratio,  and  an  indicator  flag  to  designate  either  that 
new  profile  coordinates  are  to  be  read  or  that  the  previous  section  profile  is 
to  be  repeated.  The  next  card  specifies  the  numbers  of  coordinate  pairs  in 
the  upper-  and  lower-surface  data  tables  for  the  section  and  an  indicator  flag 
to  designate  whether  the  section  is  symmetric  or  asymmetric.  Next,  if  a  new 
section  profile  is  being  given,  the  upper-surface  coordinate  pairs  of  the 
profile  are  read  proceeding  from  leading  edge  to  trailing  edge,  followed  by 
the  lower-surface  coordinate  pairs  in  the  same  order.  If  the  profile  at  the 
new  section  is  identical  to  that  of  the  previous  section,  surface-coordinate 
data  are  read  internally  by  the  program  and  need  not  be  repeated  in  the  input 
data  list.  Also,  if  a  section  profile  is  symmetric,  only  the  upper-surface 
coordinates  are  required  in  the  input  data  set;  the  lower-surface  coordinates 
are  set  within  the  program. 

The  code  provides  the  following  tabular  output: 

•  input  data  listing, 

•  computational  coordinates,  stretching  derivatives  and  grid  benchmark 
values , 

•  wing  configuration  data  at  grid  stations, 

•  maps  of  surface-adjacent  grid  points  (upper  surface,  lower  surface, 
side  edge  in  planform  view), 

•  iteration  summary, 

•  final  circulation  distribution, 

•  detailed  solution  arrays  (if  requested),  and 

•  pressure-coefficient  and  surface-slope  distributions. 

The  input  listing,  coordinate  tables,  and  wing  data  provide  a  check  of  problem 
set-up.  The  grid-point  maps  define  the  wlng/body  configuration  in  the 
computational  domain.  The  iteration  summary  printed  after  each  complete 
relaxation  cycle  lists  the  value  and  location  of  the  maximum  potential 
Increment  between  the  current  and  previous  iterations,  the  number  of 
supersonic  points  detected,  the  current  wing-root  value  of  circulation,  and 
Information  about  the  embedded  Iteration  process  required  to  update  the 
Tref f tz-plane  boundary  condition.  The  solution  arrays,  when  requested  via 
input  parameter  NPRINT,  include  the  complete  potential  array  as  well  as 
control-point  and  image-point  potential  values  which  arise  in  the  numerical 
application  of  the  surface  boundary  condition.  These  array  print-outs  are 
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useful  mainly  for  code  diagnostic  purposes.  Finally,  in  addition  to  the 
pressure-coef f icient  distribution,  both  prescribed  surface  slopes  and  velocity 
slopes  calculated  from  the  conv<^rged  potential  field  are  printed  to  provide  an 
accuracy  check  on  the  surface  boundary  conditions. 
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APPENDIX  B.  LISTING  OF  COMPUTER  PROGRAM 


P«OG«AM  TWH3(  INPUT. OUTPUT  t  TAPE baliVPUT  fTAPtfesOUTPUT) 

SOLVES  THt  FULL  POTENTIAL  EQUATION  FOR  TRANSONIC  FLOW  PAST  A 
GENERAL-WiNG/SEMl-INFINiTE-FUSELAbE  CONF IGURATION. 


C 


i 


c 


c 


DIMENSION  TITLE(20) .XTES(5) 


COMMON/CONST/ 
COMMON/DAT  AX/ 


COMMON/SOLVO/ 


ALPHA.ZMACH.OPLIM.EPSI  .wE.WG.eSPAN.  Ar«;.CSA.RXl  ,RX2. 
SNA.TOETA.TOXI.TUZETA.DETA.DXI .OZETA. ILL . IM. I  NOSE . 
ITAlL.iTE. JM« JROOT. JTIP.KM.XW. JBl » JB2. JB3.PLCB1 ♦ 

Ab.ETAO.NET  A.NXI.NZETA.XCAPO.X  10.  YCAPO. 
NF.XF (ISO) .ZF< ISO) .PF ( 150) . JSEC T . JBL ♦ J8T . NSU ( 5) ♦ 


NSL (5) .YS(S) .XLES<b) .ZL£S(5) .CS(5).ATS(S).TS(S). 


FS(5)  .XSUdSO 
-  ‘  [49.19  - 

•  ♦I 


Bbil49. 


(150.5)  .XSLdSO. 

i  lBt^AX.l5^MAX?fi 


.5) »Z5U{ 150.5) 
(49) 
amma  ( 


.ZSL (150.5) 


JMAXl* jMAXT.KMAXl*KMAXT.NSUP.P(49.19.25) .PLEl (19) . 
PNew(25) . PNOSE. PS TM (49.25) . PT 1 d 9. 25) . PT2 (1 9 . 25) . 

PT 3  (19.25)  .PTEl  (19)  .PTE2d9)  »PwbL  (49.19)  .PWBU(49«  19) 


<TITLE<N) ♦N=1.20) 

FNXI , Al .A2.XI0.XCAP0 
FNETA. A3.A4.ETAO. YCAPO 
KNZETA. AS 


MSTART=0 
NCYCLE=0 
READ  (5.901) 

READ  (5.902) 

READ  (5.902) 

READ  (5.902) 

NXlsFNXl 
NETAaFNETA 
NZETAaFNZETA 

READ  (5.902)  F  I  TERM ♦ ENH ALE .F NPR I  NT 
ITERMaFITERM 
NHALFrENHALF 
NPRINT=FNRRINT 

READ  (5.902)  4E,WG.DPLIM,EPSI 
‘  ZmaCH. ALPHA 

(TITLE (N) .Nal .20) .ZMACH, ALPHA 
ITERM.NHALF.NHRINT  .*»£,WG,DPlIM.EPSI 
NX!.Al.A2.XIO,xrARO.NETA,A3.A4.ETAO' 


READ  (5.902) 


WRITE 

WRITE 

WRITE 


(6.9031 

(6.904) 

(6.905) 


vr'Anrv.  m7Ct  a  .  au 


READ  (5.902)  FNF 
NFaENF 

IF  (NF.EQ.O)  60  TO  150 
DO  llO  Nal.NF 

no  READ  (5.902)  XF  (N)  .ZF  (N)  .RF  (N) 
WRITE  (6.911) 

DO  120  Nal.NF 

120  WRITE  (6.912)  N« XF (N) . ZF (N) .RF (N) 
150  CONTINUE 


READ  (5.902)  FJ5ECT.FJHL.FJBT 

JSECTaFJSECT 

J0L=FJ8L 

JRTaf jbT 

J  =  0 

211 


i 

! 


i 


I 
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IF  ( J.GT.  JSECT)  GO  TO  f?37 

READ  (5»y02)  rs(J)  .XLES(J) .ZLES(J> .C5(J) tATSI  J>  .TS(J) *FS(J) 

IF (FS(J) ,tU.O. )  GO  TO  230 

READ  (5»<J02)  FNSIJ.FNSL«FKSYM 

NSU ( J) =FMSU 

NSL  ( J)  =FNi>L 

ksym=fksym 

NU=NSU( J) 

DO  218  I=1»NU 

READ  (!3«902)  xsu  ( I  .  J)  ♦  ZSU  ( I  ♦  J ) 

218  CONTINUE 
NL=NSL<J), 

IF  (KSYM.EO.l)  GO  TO  220 
00  219  I=1*NL 

READ  (5.902)  XSL ( I » J) »ZSL ( I » J) 

219  CONTINUE 
GO  TO  211 

220  00  222  1=1. NL 
XSL ( I . J) =XSU ( I . J) 

2SL ( I . J) =-2SU ( I . J) 

222  CONTINUE 
GO  TO  211 
230  NSU(J)=NU 

DO  235  1=1. NU 
XSU(I.J)=XSU(1»J-1) 

ZSU(I.J)=ZSU(I.J-1) 

235  CONTINUE 
NSL(J)=NL 

00  236  1=1. NL 
XSL (I . J) =XSL ( I . J-1 ) 

ZSL(I.J)=ZSL(I.J-1) 

236  CONTINUE 
60  TO  211 

237  CONTINUE 

WRITE  (6.921)  JSeCT.JbL.J8T 
00  380  Jsl.JSeCT 

WRITE  (6.922)  J . YS < J » ♦ XLES ( J ) . ZLES ( J ) .CS ( J ) *  ATS ( J) * TS ( J) . FS ( J) 
IF  (FS(J) .NE.O.)  GO  TO  370 
WRITE  (6.923) 

GO  TO  380 
370  NUsNSU(J) 

WRITE  (6.924)  NU 

WRITE  (6.926)  ( XSU < I . J) . ZSU ( I . J ) . 1 =1 . NU ) 

NLsNSL ( J) 

WHITE  (6.925)  NL 

write  (6.926)  (XSL(I.J) .ZSL(I.J) ♦i=l*NL) 

380  continue 

00  A38  J=1.JSECT 

438  XTES (J) =XLES ( J) ♦CS ( J) 

IF  (JBL.Nt.O)  GO  TO  440 
JSECTX=JSECT-1 

XSLOP* (XLES( JSECT) -XLES{ 1 ) ) / ( YS ( JSECT ) -YS (  I ) ) 
ZSL0P=(2LES(JSECT)-ZLES(1) )/(YS(JSECT)-YS( 1) ) 

00  439  J=2.JSECTX 

XLES ( J) =XLE5( 1 ) ♦ ( YS( J)-YS( 1 ) ) *XSLOP 

439  ZLES(J)=ZLES(1) ♦ ( YS ( J ) -YS ( 1 ) )*ZSLOP 
GO  TO  445 

440  JBLXsjBL-1 

XSLOP=(XLES(JbL)-XLES(l)  )  /  ( YS  ( JttL  ) -YS  ( 1 )  ) 
ZSLOP=(ZLES(JbL)-ZLES(l) )/(Y5(JBL)-YS(1) ) 

DO  441  J=2.JBLX 

XLES ( J) =XLE5 (  ) ♦ ( YS ( J) -YS ( 1 ) ) *XSLUP 

441  ZLES ( J) =ZLES ( 1 ) ♦ ( YS ( J) -YS( I ) ) •ZSLOP 
JBLP=JBL*1 

J5ECTX=JStCT-l 

X5L0P=(XLES( JSECT) -XLES (JBL) ) /( YS ( JSECT ) -YS ( JHL ) ) 

ZSL0P= (ZLES (JSECT) -ZLES (JRL) ) /( YS (JSECT ) -YS (JBL ) ) 

00  442  JsJBLP* JSECTX 

XLES (J)=XLES (JBL) ♦ (Y5(J)-YS( JHL) ) *XSLOP 
44?  ZLES (J)=2LES( JBL) ♦(YS(J)-YS( JBL) ) *ZSLOP 
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c 


c 


445  IF  (JRT.Nt.O)  GO  TO  4*7 

XSLOP=  (XTE.S  ( JsFC  T)  -xTtS  ( 1  )  )  /  ( Yb  ( JSE.C  T  )  -YS  (  1 )  ) 

00  446  J  =  «:»JSECTX 

446  XTES(J)=XTES(1)* (YS( J) -YS( 1 ) ) *XSLUP 
GO  TO  460 

447  JHTX=jHT-l 

XSLOP= (XTtS  < JdT ) -XTFS(1))/(YS(JHT) -YS ( 1 ) ) 

DO  44a  JsiftJrifx 

44B  XTE6 ( J) =XTeS( 1 ) ♦ (YS(J> -yac 1 ) )*XSLOP 
J8TP=JtiT*l 

XSLOP=<XTt;S(JdECr)-XTES(JHT)  )  /  (  YS  (  JSECT  )  -  YS  (  JBT  )  ) 
DO  449  JaJdTPi jSeCTX 

449  XTES(J)=XTES(JBT)«(YS(J>*YS<JRT) )*X5L0P 

450  continue 

00  451  J=ltJSECT 

451  CS(J)=XT£S(J)-XLES(J) 
dSPAN  =  £*YS  < JStCT 1 


510 


520 


550 

590 

600 


901 

90? 

903 


904 


905 


4PITF  (6»931) 

ITE81=1 
GO  To  520 

wPITe  (6*932) 

IT£H1  =  1TE><M 

CALL  hale? ( I TEHM.NPRINT ) 
ITER2=ITEW1*I Tepm-1 
CALL  6«I0 
GEOM 
PHOEL 

in! T (MSTART) 

(6*941 ) 


CALL 

CALL 

CALL 

CALL 

write 


DO  590  ITEhsI  TERI *ITFR2 
CALL  SOLVt(lTER) 

IF  (NPHINT,E(J.0)  60  TO  550 

IF  (  !  T£«/NPWlfvr*NPRlNT.NE.  I  TER)  GO  TO  550 

CALL  ARRAYS (NPRINT ) 

CALL  CPCOMP 

IF  (OP'aAX.LE.UPLlM)  GO  TO  600 

continue 

CALL  ARRAYS (NPrINT ) 

CALL  CPCOMP 

IE  (NCYCLt.EQ.NHALF )  STOP 

MSTARTsI 

go  to  510 


format 

format 

FORMAT 


(20A4) 

(HE  10.5) 
(33H1TRANS0NIC 


WING/dODY  program  TwB3 


G.  E.  LHMIELEwSKI 
MCDONNELL  DOUGLAS  Ml 
ST,  LOUIS*  MISSOURI 
JUNE  1980 


format 


format 


52X43HC00E0  PYi 
85X43H 
85X43H 
85X43H 
/1X20A4/// 

5X40HMACH  number 
5X40HANGLE  OF  ATTACK 
3X7H0EGREES////) 

(21H-EXECUTI0N  parameters// 
5A40HMAXIMUM  ITERATIONS 
5X40HGRIO  HALVING  CYCLES 
5X40HSOLUT JON  ARRAY  PRINT  CYCLE 
1X21HRELAXAT ION  PARAMETERS// 
5X40HELL IPTIC  POINT 
5X40HCIHCULATION 
5X40HCONVERGENCE  LIMIT 
5X40HUAMPING  FACTOR 
(22H-CU0RDINATE  PARAMETERS// 
5X40H5TRFAMWISE  DIRECTION 
5X40H 
5A40H 
5X40H 
5A40H 


[SEARCH 

63166 

LABS./ 

/ 

/// 

2MACH 

• 

E8. 3/ 

ALPHA 

s 

E8. 3* 

ITERM 

2 

18/ 

nhalf 

s 

18/ 

NPRINT 

= 

18//// 

RE 

s 

E8.3/ 

wG 

s 

E8.3// 

<PLIM 

s 

E12.3/ 

EPS! 

s 

E8.3/) 

NXI 

• 

18/ 

A1 

3 

F8.3/ 

A2 

8 

E8.3/ 

AIO 

s 

F8.3/ 

ACAPO 

3 

E8.3// 
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or-o  nr>or>r»r> 


"1 


« 

5X40HSPANWISE  DIRECTION 

NETA 

s 

18/ 

« 

5X40H 

A3 

3 

F8.3/ 

5A40H 

A4 

5 

F  8.3/ 

5X40H 

ETAO 

S 

F8,3/ 

» 

bX40H 

YCAPO 

5 

F8.3// 

» 

5X40HVFRTICAL  DIRECTION 

NZETA 

3 

18/ 

# 

6X40H 

AS 

3 

Fe.  3) 

911 


format 

format 

format 


« 

* 

922  format 

* 

* 

* 


921 

924 

925 
92f, 
931 


format 
format 
format 
forma  t 
format 


» 

932  format 

« 

* 

format 


94] 


* 


(  1H-///1X1  3HFUSELAOF.  DA  T  A// /9X  IHN*  8X2HXF  ♦  8X2MZF  f  8X2HRF/ ) 
(I10t8F10,4) 

(lH-///inH  WIND  OATA///20H  NUMBER  OF  SECTIONS  *110/ 

22H  SECTION  AT  L.E.  RREAK»I8/ 

22H  SECTION  AT  T.E.  HRtAK»l8) 

(1H-///13H  WIND  SECTION. 1 10// 

1X5HYS  =  F9.6.5X5HXLES=  FV.B.SXSRZLES*  F9.6. 

SXShCS  =  F9.&.SXSHATS  =  F9.6»5X5Rrs  *  F9.6. 

SXSHFS  =  F9.B) 

(43H-SORFACE  OATA  IDENTICAL  TO  PREVIOUS  SECTION) 

(1H-/26M  UPPER  SURFACE  (X-Z  P A I RS ) // 1 X5HNSU  =  18/) 
<1H-/20H  lower  surface  (X-Z  PAIRS) //IXSHNSL  *  18/) 

<  UFl  1 .6) 

<1M-///4SH  INITIAL  start  *****4********4******#**##**** . 
51H**<******«**«*»*********************4******** 
40H»**«»«*«««»*«*»»»**«******«*»«»**»»*»»*«  > 

(1H-///4SH  CONTINUATION  ON  HALF-SPaCED  GRID  *•*♦***♦*♦*, 

51H**************************************************** 
4 •**«»******«*•*****•** •♦•*«*♦*****♦♦*  ) 

( lH-///lX17HlTtRATION  SUMMARY// 

6X4H1  TFR.  10X5HDPMAX.4X1MI  .4XHJ.4X1HK  .6X4HNSUP. 
14X6HGAMMAR.5X5HITERT. 7X6 10PMAXT. 3X2HJT .BXZHKT/) 


END 


subroutine  HALFSdTERM.NPHiNT) 

:IRCULAT10N  ARRAYS  ONTO  A  HALF-SPACED 


interpolates  potential  and  Cl 

COMPUTATIONAL  GRID. 


COMMON/CONST/  ALPHA.ZMACH.OPLIM.EPSI ♦WE.WG.0SPAN.Al2.CSA.RXl »RX2. 

SNA.TDETa.TOXI.TOZEU.OETa.OXI.OZETA.ILE.IM.INOSE. 
ITAIL.ITE.OM.UROOT. JTIP.KM.XW. J81 . JB2. J03.PLCB1 . 

COMMON/OATAX/  . A5.ETA0.NET A.NX I .NZETA.XCAPO.XIO.YCAPO. 

NF. XF < 1 50 )♦ ZF ( 1 50 ) .RFI 1 50) » JSECT » JBL ♦ J0T . NS J (5) . 
NSL(5) .YS(5) .XLEST5) .ZLESIS) »CS<5) .ATS(5) .T5(5) » 

FS(5)  .XSU(I50.5>  fXSLdSO.S)  .ZSyilSO.S)  .Zr  - 

COMMON/SOLVO/  OL  I  <*9.  l9)  .DL2U9.I9)  .DSl  {49)  .052  (49: 

DU2(49»19) ♦0PMAX»UPMAXT.6AMMA(I9) .JMi .  _ 

PT5d9.25)  .PTEU19)  .PTE2d9>  .PW0L(49.  I9)  «PW8U(49.l9) 
DELL <49. 19) .DELS (49) .OELU (49.19) .DZDXL (49. l9) . 
0ZDXU(49.19) .0ZDYL(49. 19) .OZOYU (49. 19) »MWaL(49.19) . 
HWBU<49.19)  .IZLd9)  .  IZU(19)  .JBOO<49)  ♦K800L(49.  19)  . 
K0ODU(49.19) .ZL(49.19) .ZU(49»19) 


. I  9 t 

2SL(150.5) 

Vlldfi”' 


common/surf  / 


NXla2.*NXl 

NFTA«2.«NETA 

NZeTA*2.*NZETA 

ITeRM«ITeRM/2 

nprint*2*nphint 

IMNa2*IM-l 

jMN=2*JM-i 

KMN»?4KM-1 
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c 


c 


c 


IMNX=IMN-i 
JMNXr JMN-1 
KMNXsKMN-1 
DO  ^00  I=lfIM 
JE  =  JPOO  (  I ) 

JtX=JE-l 
DO  1,10  JsJEfJ'^ 

DO  110  1^  =  1. KM 
KN  =  KMN- ( 2*K-1 ) ♦  1 
110  P(I»J.KM=Pn«J,K) 

DO  120  K=2tKMNX,2 

120  P(1»J»K)=0.5*(P(I»J.k-1)*P(I,j,K*1)) 
130  CONTINUE 

IF  (JE.EO.l)  DO  TO  20U 
DO  190  J=1*JEX 
KE=KRODU (I . J) -1 
P ( I » J»KE) =0U1 ( I . J) 

00  140  K=Kw,KM 
KN*KMN- (2*K-1 ) ♦! 

140  P( I » J»KN) =P ( 1 tJ.K) 

KWN1=2*KW 

DO  150  K  =  K*<N1  ♦KMNX«2 
150  P(lfJ»K)=U.5*(P{ItJ»K-l)*P(I,JtK*l)) 
KE=KROOL ( 1 » J) ♦ 1 
P( I »J.KE) =DL1 ( I *J) 

00  160  K=1*KW 
KN=KMN- <  2*K-  1  ) ♦  1 
160  P ( I » J»KN) =P ( I » JtK) 

KWNl=2*KW-2 
DO  170  K  =  2»kWim1,2 

170  P(I.J»K)=0.5*(P(NJ»K-1)*P(I,J*K*1)) 
190  CONTINUE 
200  CONTINUE 


KWNs?*K»(-1 
DO  210  I=INOStfIM 
JEXsjROD ( I) -1 
210  P(I.JEX»K»iN)=OSl  (I) 

PdNOSE*  1«KWN1=PN0SF 
DO  220  JaltJTIP 
P<lLE»J»K»iN)sPLEl  (J) 

220  P(  ITE  ,  J.K*«N)  =PTE1  ( J) 

00  260  K=1»KMW 
00  250  J=ltJM 
no  230  I=lfIM 
IN=IMN- (2*1-1 ) ♦! 

230  P ( IN, J.K) sP< I . J,K) 

DO  240  I=2,IMnX,2 

240  P(I,J,K)=0,5*(P(I-l»JtK)*P(I»l»J»K)) 
250  CONTINUE 
260  CONTINUE 


310 


320 

330 

340 


P (2*IN0SE“1 f 1 *KWN) =0S1 ( I NOSE) 

P (2*ILE»JTIp«KwN) =0S1 (ILE) 
P(2*lTE»JTIP,KWN)rOSl(ITE) 

DO  340  K=l,KMrN 
DO  330  lalfIMiN 
DO  310  J=1*JM 
JN  =  JMN- (2*JM-i ) ♦  1 
P(I,JNfK)=P(I,J,K> 

DO  320  J=2,JMNX,2 

P(I,J,K)=0,5*(P(I,J-ltK)*P(I,J«’l»K)) 

CONTINUE 

CONTINUE 


DO  410  J=1,JM 
jNxJMN- (2*J-l ) ♦  1 
410  gamma (JN) =GAMMA (J) 

DO  420  J  =  2fJMtNX,2 

4  20  gamma (J) =0.5* (GAMMA (J-1 ) ♦GAMMA (U^ I  )  ) 
JT1PN=2*JTIP-1 
gamma (JT IHN4 1 ) =0. 
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JR00TN  =  «?*JH00I-1 
gamma  (JROOTn-<!)  =0. 

return 

END 


subroutine  grid 

CALCULATES  GRID-POINT  COORDINATES*  STRETCHING  OERIVIATIVES*  AND 
WING  CONFIGURATION  DATA. 


dimension  XI (49) . YCAP ( 19) 
COMMON/CONST/ 


COMMON/COORD/  CHORD < 1 9) . OCOY ( 1 9) ♦0TDY()9) »DXLEDY( 

.  e-t  A  kAi  c  .  ^  a  kA  hi  ^  fl  a  kAC  ^ 


♦FLE*FN.Ft 
TaU(19) .X( 


COMMON/DATAX/ 


(19) ♦0TDY()9) »DXLEDY(19) ,ETA(1?) .F(4?) 
19) ♦6AMLE*GAMN.6AMF*6AMTE*GAMWt»H(251 ♦ 
9)  ♦XCAP(49)  ♦XLE(19)  fYIW)  ♦2ETA(25)  »6F. 

f\*9i  .  Ai  r4LJAV#%A«  r\*r\w/t04  .AaijT 


lAUli^  /  fT  \  fvr  i 

ZCAP(25) .ZLE(I9) ♦0ZLEDY(19) * ALPHAT ( 19) ,DADY ( 19) tGwT 
Al «A?«A3«A4«AS*ETa5«NETA«NXI«N2ETA«XCAPO»XIO*YCAPO* 
NF*XF(150)  ♦2F(150)  tRFdbO)  ♦  JSECT  *  JBL*  J8T*NSU  (5)  « 

NSL (5) ♦YS(5) »XLES(5) .ZLES(5) *CS(5) *415(5) *TS(5) * 
FS(5) *XSU(150*5) «XSL(150*5) *ZSU(150*5) *2SL(150f5) 


102  format 
1 
2 

108  format 

294  format 
1 
2 

355  format 

357  format 
1 

411  format 
2 

595  format 

655^E0RMaT 

1 

657  format 

734  format 
1 
2 

7749  format 

BOl  format 
036  format 


3 

4 

5 

6 
7 

871  format 
1 
2 


( 1M-///1X27MSTREAMWISE  GRID  COORDINATES/// 

1X8HDXI  s  F14,6// 

9X1HI *  1 1X4HXCAP* 13X2HXI * 14X1HF/) 

(I10*3(1X*F14,6) ) 

(IHO/IXHHFLE  =  F14.6/1X8H6AMLE  *  F14.6/ 

1X8HXCAPLE  *  F14,6/IX8HX1LE  =  F i4.b//lX8HFTE  »  F14.6/ 
1X8HGAMTE  =  F14,6/IX8HXCAPTE  s  F 1 4 , 6/ 1 X8HX I TE  *  F14.6) 
(//IX8HFM  s  F14.6/1X0HGAMN  s  F  1 4 . 6/ 1 XSHXC APN  »  F14,6 

/1X8HXIN  s  F14,6) 

(//1X8MINOSE  =  I14/IX8HILE  =  Ii4/1X8HITE  »  114/ 

1X8HITAIL  =  114/lXeHIM  =  114) 

( 1H-///1X25HSPANWISE  GRID  COORDINATES/// 

1X8H0ETA  a  F14.6// 

9X1HJ* 1 1X4HYCAP* 12X3HETA* 14X1MG/) 

(1H0/1X8HGWT  =  F14.6/1X8HGAMWT  »  F14.6/ 

1X8HYCAPWT  =  F14,6/lX8MtTAWT  =  F14.6) 

(//1X8MGF  a  F14,6/IX»HGAMF  a  F 1 4 . 6/ 1 X0HYC APF  a  Fl4,6 

/IXSHETAF  a  F14,6) 

(//1X8HJR00T  a  I14/IX8HJTIR  =  I14/1X8HJM  »  114) 

(1M-///1X25HVERTICAL  GRIU  COORDINATES/// 

1X8H02ETA  a  F14,6// 

9XIHK* 11X4HZCAP*! 1X4HZETA* 14XIHH/) 

(//1X8HKW  a  I14/1X8HKM  r  114) 

(I8*10(1X,F10.6) ) 

( 1H-///1X35HWING  PLANFORM  AREA  S  *  F8.4/ 

1X35HWING  ASPECT  RATIO  AR  a  F8.4/// 

IX35HSTATI0N  AT  L.E,  BREAK  YBl  *  F8.4* 

10X5HJB1  X  I4.i0X7HETABl  *  F8 , 4 *  I UX 7HPLCR 1  *  F8.4/ 
1X35HSTATION  AT  T.E.  BREAK  YB2  a  F8.4* 

1UX5HJH2  a  I4,10X7HeTAd2  a  F8 , 4 *  1 UX 7HPLCB2  »  T8.4/ 
iX35HCONSTANT-CHORO  STATION  YB3  a  F8.4* 

10X5HJB3  a  1 4 1 1 0 X7HE T AB3  =  F8 , 4 *  1 0X7HPLCS3  *  F8.4) 
(1H-///1X41HWING  CONFIGURATION  DATA  (PHYSICAL  DOMAIN)// 
7X1hJ*8X3HXLE*4X7HDXLE/UY*6X5HCH0RD*6X5H0C/DY*8X3HTAU* 
6XbH0r/0Y,8X3MZLE*4X7HD2Lt/DY»5X6HALPHAT*6X5HDA/0Y/) 
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oo 


C 


XCAP-STRE  rCHiNG 

HPI=2.*ATAN ( 1 . ) 

X IM=1 . ♦XIO 
0X1=2. •XIM/NXI 
IM=NXI*1 
I0=NXI/2*1 
I0=XI0/DX1 
lLt=I0-I0 
ITE=I0*I0 
XI ( I0)=0. 

I  I  =  I04l 
00  50  I  = 

50  XI  (  I ) =XI ( I-l ) *0X1 
DO  t)0  I=«?»I0 
IX=I0-I*1 

50  XI (IX)=XI (IX^I)-OXI 

ALCU0.5*  (3,*XCAP0/X10-HPI*A1  ) 

ALC2=0.5*  <HPI»a1-XCAPO/XIO)/ (XI0*XI0) 

WRITE  (6»102»  DXI 
XCAPIDs-l.E^^Q 
F ( 1 ) =0. 

ILX=ILE-1 

00  110  I=1«ILa 

IF  (I.EQ.l)  GO  TO  109 

PP  =  XI  { I ) ♦XIO 

TiM2  =  TAN  (HPI«PP) 

PP3=PP**3 

TfvJ3  =  TAN  (HPI*PP3) 

XCAP (I ) s-XCAP0^Al*TN2^A2*TN3 

F ( I ) =1./ (HPI* (Al* ( I.4TN2*TN2) ♦3.*A2*PP*PP* ( 1 , ♦TN3*TN3) ) ) 
109  WRITE  (6*100)  I tXCAP ( 1 ) ,XI ( I ) .F ( 1 ) 
no  CONTINUE 

DO  120  I=ILE«ITE 
PP2=XI ( I ) 

XCAP(  I )  sXKl )  *  (  ALC1*ALC2*PP2) 

F(I)=l./(ALCn3.*ALC2*PP2) 
write  <6fios)  i*xcAP<i) »xi ( n .F ( n 
120  CONTINUE 

XCAP (IM) =1.E^99 
F(IM)aO. 

itep=ite*i 
00  130  i=Itep*im 
IF  (I.EU.IM)  GO  TO  129 
PP  =  XI  ( I ) -XIO 
TN2=TAN(HPI*PP) 

PP3=PP**3 
rN3=TAN (HP1*PP3) 

XCAP ( I ) =XCAP0*A1*TN2^A2*TN3 

F(I)=1./(HPI*(a1*( 1 .♦TN2*TN2) ♦3.*A2*PP*PP* ( I .♦TN3*Tn3)  )  ) 

129  WRITE  (6*100)  I  * XC AP ( I ) , X I ( I ) » F ( I ) 

130  CONTINUE 

IF  (XCAPO.NE.0,5)  GO  TO  200 
xile=-xio 

XITE=XIO 
GO  TO  290 
200  XI1=XI0^UXI 
210  PP  =  xn-XlU 

TN2=TAN (HPI«PP) 

PP3=PP**3 
rN3=TAW (HPI*PP3) 

FUN=A1*TN2>A2*TN3-0.5^XCAPO 

FUNPH=rtPl* (Al* ( 1 .♦TN2*Tn2) ♦3.*A2*PP*PP* ( I .♦Tn3*TN3) ) 

OEL=FUN/FUNPR 

Xl2  =  Xll-(JtL 

IF  (OEL.LT, 1 ,E-10)  GO  To  250 
XI1=XI2 
GO  TO  210 
250  XITE=XI2 
xtLE=-XI2 
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nno 


C 


29n  <iAMLF  =  ABS(XILfe.-XI  ( ILF-1)  ) 

GAMTE  =  ABS(XI ( I TE ♦ 1 ) -X I  TE ) 

PP=XIT£-XI0 
TN2=TAN (MPI*PP) 

PP3=PP**3 

TN3  =  TA;vI  (HHI*PP3) 

XCAPTE=XCAP0»AI*TN?*A2*TN3 

XCAPLE=-XCAFTE 

FTE=1 ./ (HPI* ( Al* ( 1 .♦TN2*TN2) ♦3.*A2*PP*PP* ( 1 .♦TN3*TN3) ) ) 
FLE=FTE 

*)PITE  (6.294)  FLF.6AMLF.XCAPLEtXlLE.FTE.GAMTE.XCAPTE.XITE 


DCS=CS(2)-CS(1) 
l)XS  =  XLEb(2)-XLES(l) 

DYS=YS(2)-YS( 1) 

IF (DYS.EO.O, )  STOP  tFPPOR  YS ( 1 ) = TS ( 2) • 

XREF  =  XLES(  1 )  -  YSd  )  •UXS/DYS 
OXS=OXS*OCS 

CREF=XLES( 1 ) tCS ( 1 ) -YS ( 1 ) *DXS/D YS-XREF 
XCAPN=-XRtF/CREF-0.5 
no  300  1  =  1  . ILE 
OIFF=XCAP(I)-XCAPN 
IF  (OIFF.GE.O.)  GO  TO  305 
300  CONTINUE 
305  INOSE=I 

XI1=XI (INOSE-1) 

310  PPsxntXIU 

TN2=TAN (HPI*PP) 

PP3=PP**3 
TN3=TAN (HPI*PP3) 

FUN  =  Al*TN2'tA2*TN3-XCAPN-XCAP0 

FUNPRsMPI* (Al» ( U*TN2*TN2) *3. •A2*PP*PP* ( 1 .♦TN3»TN3) ) 

OELsFUN/FONPR 

XI2=XI1-DEL 

ueLTsABS(OEL) 

IF  (OELT,lT,1.E-10)  go  to  350 
XI1=XI2 
GO  TO  310 
35o  X1N=XI2 

GAMNsABS  (XI  (  IinOSF-1  )  -XIN) 

PP=XIN>XIO 

TN2=TAN(HPI*PP) 

PP3=PP**3 

TN3=TAN(HPI*PP3) 

FN=1 ./(HPl* (Ai*( I .♦TN2*TN2) ♦S. *A2*PP*PP* ( 1 . ♦TN3*TN3) ) ) 
WRITE  (6.355)  FN.GAMN.XCAPN. XIN 
FL  =  XF (NF) -XF (  I ) 

XCAPTx(FL-XREF) /CREF-,5 
DO  360  1  =  1. IM 

IF (XCAP( I) .GT.XCAPT)  GO  TO  370 
360  CONTINUE 
370  ITA1L=I-1 

IF(XF(NF) .GT,25.)  ITAIL=IM 

WRITE  (6.357)  INOSE . ILE . I TE ♦ I T AIL .  IM 


ycap-stretching 


ETAM=l .tETAO 
UETA=ETAM/NETA 
JMxNETA. 1 
JTIP=1.^ETA0/UETA 
ETA  (  1 ) =0. 

DO  400  J=2.JM 
400  ETA(J)=ETA(J-i) tOETA 

ALC3=0.5* (3,*YCAPO/ETaO-HPI*A3) 

ALC4a0.5* (HPI«A3-YCAP0/ETA0) /(ETA0«ETA0) 
WRITE  (6.411)  OETA 
00  420  J=1.UT1P 
PP2=ETA ( J) **2 

YCAP( J) =ETA(J)* ( ALC3tALC4*PP2) 
G(J)=1./(ALC3^3.*ALC4*PP2) 

WRITE  (6.108)  J.YCAP(J)  .ETAU)  .G(J) 
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4^n 


44^ 

450 


500 

510 


550 

590 


C 


460 

600 

605 

610 


650 


C 

C 

C 


CONTINUE 

YCAH  (  JM)  =1  .E^9‘) 


(>(  =0. 

DO  450  J=JTIPP,JM 
IF  (J.EQ.JM)  00  TO  445 
PPafcTA(J)-ETAO 
TN^sTAN(HPI•PP) 

PP3=PP**3 

TN3=TAN(HPI*PP3) 

YCAP( J) =YCAP0*A3*TN2*A4*TN3 

Q( J) *1 ./ (hPl«  U3« ( l.*TN2*TN2) ♦ J . * A4*PP#PP* < 1 , ♦ TN3* TN3) ) ) 
WPITE  (6*108)  J*YCAP ( J> «ETA ( J) *6( J) 

CONTINUE 


IF  (YCAPO.NE.0,5)  GO  TO  500 

ETAWT=£TAU 

GO  TO  590 

ETA1=ETA0^D£TA 

PP=ETAi-ETAO 

TN2=TAN(HPI*PP) 

PP3=PP**3 

TN3aTAN(HPI*PP3) 

FUNsa3*TN2^A4’*TN3-0.5*YCAP0 

FUNPpsHPI* (A3* ( 1 .♦TN2*TN?) ♦a, *A4*PP*PP* ( 1 

OELsFUN/FUNPR 

ETA2=ETA1-UEL 

IF  (OEL.LI . l.E-10)  60  TO  550 

ETA1=£TA2 

GO  TO  510 

ETAWTsETA2 

GAMwT=ABS(ETA(JTIP^1)-ETAWT) 

PPsETA»lT-tTAO 

rN2®TAN(HPI*PP) 

pp3spp**3 


.♦TN3*TN3) ) 


TN3STAN (HPI*PP3) 

YCAPWTsYCAP0*A3*TN2>A4*TN3 

GWTsi ,/ (HPI* (43* (1,»Tn2»TN2) ♦3,*A4*PP*PP* ( l*^TN3*TN3) ) ) 
WRITE  (6«395)  GwT,GAMwT*yCAPWT*fcTAW7 


YFSsRF (1) 

00  460  N*2»NF 

IF (YFS.LT.RF (N) )  YFSsRF(N) 

CONTINUE 

yCAPF=YFS/bSPAN 

00  600  J=l»JTiP 

OIFFsYCAP' J)-YCAPF 

IF  (OIFF.GE.O.)  GO  TO  605 

CONTINUE 

JROOTsJ 

ETAl=ETA(JWOOr-l) 

eTAlSQ*ETAl*ETAI 

FUNsFTAl*(ALCi*ALC4*F.TAlStJ)-YCAPF 

FUNPRsalC3>3,*ALC4*FTA1S0 

0EL=FUN/FUNPR 

ETA2sETAl-OEL 

IF  (OEL.LT, 1 ,e-lO)  GO  TO  650 

ETAi=ETA2 

GO  TO  610 

ETAFsETAH 

GAMF=ABS(ETA(UR00T)-FTAF) 

GFsl ./ (ALCJ*3.*ALC4*F  TAF*ETAF) 
WRITF  (6»655)  6F .GAMF, YCAPFtETAF 
WRITE  (6*657)  JROOT ♦ JT IP* JM 

zcap-stretching 

ZETAMsI ,0 

0ZETA  =  2.*<'E7AN/NZETA 
KMsNZETA*! 

KW=NZETA/2*1 
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ZETA (KW) =0, 

KK=KW* 1 

no  710  K  =  ^K♦K■'^ 

ZEfA  (K)  =ZtTA  ((<-1 )  ♦OZETA 
00  120  K=2»KW 
KX=KW-K^ 1 

ZETA (KX) sZETA (KX*1) -OZETA 
WRITE  (6.734)  OZETA 
ZCAP ( 1 ) =-l .E*99 
H(1)=0. 

ZCAP(KM)s-ZCAP(l) 

H (KM) sO. 

00  7S0  Ksi.KM 
ZET=ZETA (K) 

IF  (K.EQ,  1 .OR.K.EQ.KM)  60  TO  749 
tMl=TAN(HPI*ZET) 

ZCAP(K) =Ab*TNl 

H(K)=l,/(Ab*MPI*(l.>Twl*TNl)  ) 

WHITE  (6,108)  K.ZCAPCK) ,ZeT,H(K) 

CONTINUE 

WHITE  (6,T749)  KW.KM 

WING  configuration  DATA 

JSL=J8L 

JSTsJHT 

IF(JSL.EQ.O)  USLsJSECT 
IF(JST.EQ.O)  JSTajSECT 
XTIP=XLES(JSECT) 

CTIP=CS(JSECT) 

SWlNGsHSPAN* (XTIP^CTIP-XREF) 

SWING=SWING-(ALES(JSL)-XREF>*YS(JSL)  ,  ,  . 

SWiNGaSWiNG- (XlIs ( JSL) -XREF*XTIP-XREF) * ( YS ( JSECT) -YS ( JSL) ) 
SWINGaSwING-(AL£S(UST) ♦CS(JST)-XRrF-CHEF)*YS(JST) 

SWiNGaSwING- (XTIP*CTIP-XLES( JST) -CS (J5T) ) • ( YS ( JSECT ) ♦ YS ( JST ) ) 
ARr8SPAN*bSP AN/SWING 

00  800  Jsl.JM 
Y(J)aBSPAN*YCAP(J) 

CONTINUE 

CALClaYS(JSECT)-YS(JSECT-l) 

CALC2aXLES (JSECT )-XLES(J5ECT-l) 

CALC3aCS (JSECT ) -CS ( JStCT-1 ) 

IF  (CALC3.E(J.0,  )  GO  TO  803 

YSET=YS( JSECT )-CS( JSECT )*CALCl/CALC3 

Y83a(YSET^YS(JSECT) )/2  ^ 

XLEB3aXLFS(JSECT) ♦ ( Y83-YS ( JSECT ) ) *CALC2/CALC 1 

XTEe3=XLES( JSECT) ♦CS (JSECT) ♦ ( C ALC3^CALC2 ) * ( YB3-Y5 ( JSEC T ) ) /CALCl 
00  802  Jal.JM 
CALCaY03-Y (J) 

IF(CALC,LT,0.)GU  to  805 

CONTINUE 

JB3=JM 

Y83=l .E*99 

PLC83»0,0 

ETA83=ETA( JM) 

GO  TO  027 
JB3=J-1 

YCAPB3=Y8J/RSPan 
ETAlaETA ( JB3) ♦UETA 
PP=ETA1-ETA0 
TN2aTAN(HPI*PP) 

PP3=PP**3 

TN3=TAN(HPI«PP3) 

FUNaA3*TN2^A4*TN3-YCAPB3^YCAP0  ^ 

FUNPRaHPI*  (A3*  ( 1  .♦TN?*TN2)  ♦3.*A4*PP*PP*  ( 1  .  ♦•TNJ*TN3)  ) 

OELsFUN/FUNPH 

ETA2rETAl-0EL 

IF  (OEL.LT.l.E-lO)  (30  TO  8826 

ETAlaETA2 

GO  TO  8025 
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Ha34 


835 

8835 


0836 


8838 


ETAB3=ET&2 

PLCB3= (ETAH3-tTA ( JP3) ) /OETA 
CONTINUE 


8039 


TO  8836 


00  828  J=1»JM 

IF  ( JBL.EQ.O)  GO  TO  830 

CALCsYb ( JbL) -T ( J) 

IF  (CALC.LT .0. J  GO  TO  831 

CONTINUE 

JBUO 

VBlsYSd) 

GO  TO  832 

YBIsYS ( J8L) 

CONTINUE 

YCAPBl=YBi/BSPAN 
IF  (JBl.Eu.O)  JSTaRT=JR00T-1 
IF  (JBl.Nt.O)  JSTA«T=JBI 
ETA1=ETA ( JSTAkT) 

ETA1SQ=ETAI*ETA1 

FUNsFTAI* (ALC3*ALCA«ETAlSQ)-YCAPBi 

FUNPRaALC3*3.*ALC4*ETAlSQ 

DEL=FUN/FUNPR 

ETA2=ETA1-0EL 

IF  (OEL.LT.l.t-lO)  GO  TO  8834 

ETA1=ETA2 

GO  TO  8833 

ETABI=ETA2 

PLCB1= (ETABI-ETA < JSTAHTI )/0ETA 
00  835 

IF(JBT,EQ.O)  oO  TO  8835 
CALCsYS(JBT)-Y (J) 

IF (CALC.LT.O. )  GO  TO  8836 

CONTINUE 

JB2  =  0 

YB2sYS ( 1 ) 

GO  TO  837 
JB2aJ-l 

YB2=YS(UBT) 

CONTINUE 

YCAPB2=YH2/BSPAN 
IF  (JB2.EO,0)  JSTaRT=JR00T-1 
IF  (JH2,NE.O)  JSTART=Jfl2 
ETAIsETa  USTART) 

ETA1S(5  =  ETAI*E  fAl 

FUNsFTAl* (ALC3*ALC4*ErAlS0)-YCAP82 

FUNPR=ALC3*3.*ALC4*ETAISQ 

OEL*FyN/FUNPR 

ETA2=ETA1-0EL 

IF  lOEL.LT , 1 ,E-10)  GO  TO  8839 
ETA1=ETA2 
GO  TO  0838 
ETABPsET  A2 

PLCB2=  <ETAB2-ETA ( JSTAHT) ) /OETa 
WRITE (8*836)  SW ING * AR ♦ YH I ♦ JB I tET AB 1 fPLCB 1 » 
I  Y82. JB2»ET AR2»PLC82»YB3» Jfi3*ETAB3fPLCB3 

SL0PFs(YS(2)-YS(l) )/(XLES(2)-<LEb(l) ) 

00  845 

XLE  (J)=XLES(1) ♦ (Y(J)-YS(l) ) /SLOPE 

OXLEDY (J)=1/SlOPE 

CONTINUE 

CALC1=YS( JSeCT ) -YS( JSECT-1) ^ 

CALC2=XLES(JSECT)-XLES(JSFCT-1) 

SL0PE1=CALCI/CALC2 

J8lP=JBl*i 

DO  847  J=JH1P*JM 

XLE (J)=XLES( JSECT) ♦ (Y(J)-YS( JSECT) ) /SLOPE  I 

OXLEOY ( J) =1/SL0PE1 

CONTINUE 
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Bb?  00  Bbi  JS=2«JbECT 

IF (YS( JS) .6E .Y ( J) )  00  TO  857 
851  CONTINUE 

DYS=YS(JSECT)-YS<JSFCT-1> 

0ZLESaZLES(J5ECT)-ZLeS(JSECT-l> 

CALC=nZL£b/OYs 

ZLEBl=ZLE5(JStCT) ♦<YR3-YS(JSECT) )*CALC 
IF (Y ( J) .GT .YRJ)  GO  TO  054 
ZEE ( J) =ZLtb3> ( Y  « J)-YB3) *CALC 
OiLEDY{J) =CALC 
GO  TO  858 
854  ZLE(J)sZLER3 
OZLEDY ( J) =0 


GO  TO  858 

857  CALCs<ZLE5(  JS>-ZLES(JS-n  )  /  (  YS  <US> -YS  /  JS-I  >  > 
ZLE ( J) =ZLtS ( JS) ♦TY ( J)-YS(JS) ) *CALC 
OZLEOY ( J) =CALC 


85fl  CONTINUE 


CALC= (CS (Z) -C5 (1 ) ) / (YS(2)-Y5( 1 ) ) 

00  888 

rF(YS<  n  *l.E.  Y<J)  )  60  TO  86Z 
CHORD!  J)  sCS  ( 1 )  •»  <Y(  J)-YS(l)  )  •CALC 
OCDY ( J) sCALC 
GO  TO  8b« 

86?  00  863  JS=Z»JSE.CT 

IF(YS(JS) .GE.Y (J) )  GO  TO  867 

863  CONTINUE 

0YS=YS<JSECT) -YS( JSECT-l) 

OCS=CS(  JSECTI-CS(vJSECT-U 

CALCaOCS/UYS 

CH0HnR3sX iEB3-XLER3 

IF (Y ( J) ,GT .Y83)  go  to  864 

CHORD (J) aCH0RUB3* (Y( J)-Y83)*CALC 

OCOY(J)sCALC 

GO  TO  868 

864  CHORD »J) aCHORUHl 
OCOY { J) =0 

GO  TO  860 

867  CALCa(CS<US)-CS<  JS-n  » /( YS  ( JS) -YS  (US-1 )  ) 
CHORO(U)sCS(US) ♦<Y(J)-YS(JS> )*CALC 

OCOY ( JJ  =CALC 

868  CONTINUE 


TROOTaTb (1) 

DO  870  JalfJSECT 
870  rS(J) s7S<U? /THOOr 

CALCa(TS(Z)-TS(lJ  J / ( YS ( 2 ) - YS ( 1 ) ) 

00  870 

IF(YS(il .LE.Y(J))  60  TO  872 
tAU(J)=TSU)  *(7  (  J)-YSU)  )*CALC 
OTDY(J)aCALC 
GO  TO  878 

87?  DO  873  JS=»2.JSECT 

IF ( YS ( JS) .GE. Y ( J) )  GO  TO  877 
873  CONTINUE 

OYS=YS(UStCT) -YS( JSECT-l) 

OTS=TS( JSECT) -TS(jStCT-l) 
CALCaOTS/UYS 

TAU83=TS( JSECT) ♦ (YR3-YS{JSECT) ) *CALC 
IF ( Y ( J) ,GT .Yfl3)  GO  TO  874 
TAU( J)=TAUB3*  <Y ( J) -YB3) *CALC 
DTOY (J) aCALC 
GO  TO  878 
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ft  TAU(J)=TAUHJ 
DTDV ( J) =0 
GO  TO  rt7H 

7  CALC= ( fS ( JS) -  I S ( JS-1 ) ) / ( YS( JS) -YS ( JS-1 I  ) 

T  AU Ij) =T  S I JS) ♦ ( Y ( J) -Yb ( JS) »  *CALG 
□  TOY  ( J) =CALC 
^  CONTTNUt 

CALC=  <ATS<2)-ATS<l))/<YS(a)-YS(l)) 

DO  88B  J=ltJM 

IF  (  Y5(  1 )  .LE.Y  U)  )  GO  TO  882 
alpha  T  (J)=ATS(l)>(Y(J)-YS(n)  *CALC 
OAOY  < J) =CALC 
GO  TO  888 

?  DO  883  JS=atJb£CT 

IF (YS(Jb) .GF. Y (J) )  GO  TO  887 

B  continue 

OYSsys ( JSECT ) -YS ( JSFC  r-1 ) 

OATS=ATb (USECT ) - ATS ( JSECT- 1 ) 

CALC=OATS/OYS 

ALPHT83=ATb (JSECT ) ♦ ( YB3-YS ( JSEC T ) ) *CALC 

IF<Y(J) .GT,yB3)  go  to  884 

ALPHA  T  (  J)  =ALPr'TR3>  (  Y  (  J)  -Yh3)  *CALC 

OAOY ( J) =CALC 

GO  TO  888 

>  ALPHAT ( J) =ALPhTH3 
OAOY ( J) =0 
GO  TO  888 

'  CALC= ( AIS  < JS) -ATS ( jS-i ) ) / (YS ( JS) -YS ( JS-1 )  ) 

ALPHA  r (J) =ATS (JS) ♦( Y (J) -YS(JS> ) *CALC 
OAOY ( J) =CALr 
I  CONTINUE 

WRITE  (8*871) 

00  898  J=i«JM 

,  WRITE  (6*801)  J,XLE(J) *OXLFDY(J>  *  CHORD (J) *OCOY(J) *TaU(J) *OTOY(J) 
1  *ZLF  (J)  *UZlEiJY  (J)  ,  ALPHAT  (J)  *DAUY  (J) 

00  900  I=1*IM 
00  900  Jsl*JM 

X  ( I  »J)  s(XCAP  ( I )  *.0.5)  *CHORn(J)  ♦ALE  ( J) 

CONTINUE 

return 

END 


SUBROUTINE  GEUM 

COMPUTES  WING/0OOY  COORDINATES  Z<X*Y)  AND  SURFACE  SLOPES  D2/DX  AND 
OZ/OY  AT  GRID-POINT  STATIONS  IN  THE  PHYSICAL  DOMAIN. 

DIMENSION  SPjlSOl.XSETdSO)  ♦S0(I50)  *RFS  ( ISO)  tZFS  ( 150)  .XFSPdSO)  t 

•  RFSPdsO)  *ZFSP<  150)  *XFSPP  ( 150)  *RFSPP(IS0)  ♦ZFSPP(i50) 

COMMON/CONST/  ALPHA.ZMACH*OPLIM.EPSI*»*E*WG*BSPAN,Al2fCSA*RXl,RX2* 

*  SNA*TDETA*TDXI.TDZETa*OETA«DXI*DZETA* ILEtIM*IN6sE* 

!  ITAIL*iTE. jM.JRO0T*JTlP*KM.KW*JBl*J82* JB3*PLCBI* 

mm  1^  I  ^  T  1  ^  n  ^ 


COMMON/COOrd/ 


CFTORO(  lbf?DCDY(l?)  *DTDY(19)  »0XLE0Y  ( 19)  tETA  ( 19)  .F<49) 
*FLE*FN*FtE*G( 19> ♦ GAMLE ♦ SAMN ♦GAMF * GAMTE tOAMWT » H ( 2b) * 
TaU( 19) *X(49* 19) *XCAP(49) *XLE ( 19) * Y ( 19) . ZET A ( 25) iGF . 
ZCAP  (25)  *Z(,E(i9)  ♦DZlEDY  (19)  » ALPHA!  <  19)  .DAOY  ( 19)  *§WT 
Ai*A2*A3.A5,A6,ETA0,NETA,NXl«NZETA,XCAP0,xi0,YCAP0* 
NF*XF(150) .ZF(I50) ♦RF(i50) « JSECT * J0L * J8T ♦ NSJ (5) . 

NSL (5) *YS(5) .XLES(5) *ZLES(5) «CS (5) *  ATS (5) tTS (5) * 
FsTB) .XSUllSOiS) *XSL(i|0*5) »ZSU( 150*5) *ZSL( 150*5) 
PELL(49,19) ,DELS(4?T,0ELU(49»19) *DZ0XL (49.  19)  . 


COMMO'4/DAT  Ax/ 


COMMON/SURF 


oloxut45!?9)  *DiDYL(49SfbV*Dz6|u(iS5?9|riHwSuU.19) 

HwBU(49.19).IZL(19) *120(19) *J0OO(49) .KBOOL (49* I9)  . 


Kfl00U(49.l9) fZL <49* 19) *ZU( 49. 19) 
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non  non 


i 

FOWMaT  ( lM-///lX?qHPLANFOPM-ADJAC£NT  GRID  POINTS 
1  //3XlHl,flx7HJaOO(I>/) 

25l  format  (IA,Ilb) 


00  50  I=1»IM 
JB00(I)=1 
00  50  J=lfJM 
^UCIf J)=0. 

ZL(I.J)=0. 

0ZDXua*J>=0, 

OZl)XL(I»J)=0. 

azDYu(i»j)=o. 

50  OZDYL ( I ♦ J) =0. 

FUSELAGE,  co-ordinates 

00  56  I=IN0SE« IT4IL 

K=I-IN0SE*1 

00  54  J=1*JR00T 

DO  5?  I l=iN0SE. ITAIL 

KlsIl-lNOae^l 

XS£T(K1)=X(I1»J) 

52  continue 

N0=ITAiL-lN0S£>l 

IFIXSEr (NO) ,Gr,XF(NF) )  N0=N0-1 

CALL  SSPLINE (NF,XF.ZF*NO»XSET»ZFS.ZFSPtZFSPP»A) 

CALL  ASPLiNE (NF  *XF  fRF,SF»NO*XStT»S0*RFS*XFSP*HFSP*XFSPP»RFSPP»2) 

IF  (K.GT.i'lO)  f^sNO 

IF(Y(J» .G£.«FS(K) )  GO  TO  55 

CALCsSORT (RFS (K) *RFS (K) -y  < J)*Y ( J) ) 

ZU( I »J) sZFS(K) ♦CALC 
ZL(ItJ)=ZFS(K)-CALC 
ORFOXFaRFSP (K) /XFSP<K) 

0ZF0XF=ZFSP(K) 

DZDXU ( I ♦ J) =DZFDXF^RFS(K) *0RF0XF/CALC 
OZOXL ( I  «J) =OZFOXF-RFS(K)*ORFOXF/CALC 
DZDYU ( I » J»  =-Y ( J) /CALC 
OZOYL  < 1  * J) =Y { J) /CALC 

54  CONTINUE 

55  J0OD(I)=J 

56  CONTINUE 

00  124  IslLEtlTE 
124  J0OD(I)sjTIp4l 

WING  coordinates 

00  150  US=1«JSECT 
NU=NSU ( US) 

DO  140  N=i»NU 
XF (N) sXSUTNf US) 

RF (N) sZSUTN, US) 

140  CONTINUE  . 

00  144  1=1LE»1TE 

XSET ( I )  =  (X ( I ♦ 1 » -XLE  < 1 ) > /CHOROT 1 ) 

144  CONTINUE 

N0=ITE-ILE*1  , 

00  146  1=ILE»ITE 
XSET(I-ILE*1I=XSFT(U 
146  CONTINUE 

CALL  ASPLlNE (NSU( JS» . XF tPF . SF tNO » XSET ♦ SO . RFS ♦ XFSP » RFSP ♦ XFSPP » RFSPP 
1  ♦Z) 

DO  14R  I=1*N0 

ZSU( ILE^I-1 tJS)=RFS(I) 

XSU ( ILE^I-1 ♦ JS) sMFSP  T I )/XFSP( I ) 
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c 


c 


14R  continue 
I5n  continue 

00  154  I=iLE»lTE 
JSECTX=JSECT-1 
no  154  JS=1»J5ECTX 
00  153  J=JP00T,JTIP 

IF(YS(l).t(J.Y(J).OR.(YS(JS).LT.r(J)  .AND.  Y  (  J)  .LE  .  YS  (  JS»  1 )  )  ) 
1  GO  TO  15^ 

GO  TO  153 

15?  CALC= (Z5U  < 1  * Jb) -ZSU( I .J54l ) ) / ( Yb ( JS) -YS ( US* 1 ) ) 

ZU  (I  .  J) =ZbU ( I » US) ♦CALC* { Y ( J) -YS ( US) ) 

ZU ( I « J) =ZU ( I » J) *CHORO ( J) 

CALC= ( XSU  U . US) -XSU ( I . JS*1 > ) / ( YS ( US) -Y5 { JS*1 ) ) 

OZOXU ( I t J) =XSU ( 1 , US) ♦ ( Y ( J) -YS ( US) ) ♦CALC 

153  CONTINUE 

154  CONTINUE 


155 


156 


157 


DO  160  JS=1,USECT 
NL=NSL (US) 

DO  155  N=ltNL 
XF (N>  =X5L IN. US) 

RF (N) =ZSL (N,U5) 

CON  r INUE 

DO  156  1=1LE.ITF 

X5ET (  I  )  =  (X (1 . 1 ) -xlE ( 1 ) ) /CHORD ( 1 ) 

CONTINUE 

no=ite-ile*i 
DO  157  1=1LE*1TE 
X5ET ( I-ILt*l ) =XSET (I) 
continue 

CALL  ASPLINE (nSL (US) . XF . RF . SF . NO ♦ XSET ♦ SO . RFS ♦ XFSP ♦ RFSP » XFSPR . RFSPP 
1  .3) 


DO  158  l=liNO 

ZSL ( ILE*I-1 .US) =RFS( I ) 

XSL  (ILE*I-1  .U3)  =«,-SP  (  I )  /XFSP(  I  > 

15R  CONTINUE 
160  CONTINUE 

DO  164  I=lL£»iTE 
JSECTX=USECT-1 
DO  164  US=1.USECTX 
DO  163  UsUWOOT.UTIP 

IF(YS(l).tQ.Y(U).OR.(YS(US).LT.Y(U) .AND. Y (U) . LE . YS  ( US*  1 )  )  ) 
1  GO  TO  16^; 

GO  TO  163 

16?  CALC= (ZSL ( I . US) -ZSL ( I *US*1 > ) / ( YS (US) -YS(US*1 ) ) 

ZL  (  I  .U) =ZSL ( I »US) *CALC*(Y( J)-YS(Ub) ) 

ZL ( 1 «U) =ZL ( I .J) ♦CHORD (U) 

CALC= (XSL ( I .US) -XSL ( I .US* 1) ) / ( YS (US) -Y5 (US*1 ) ) 

DZOXL ( I .U) =XSL ( I .US) * ( Y (U) -YS (US) ) *CALC 

163  continue 

164  CONTINUE 


DO  17?  1=ILE.1TE 
on  17?  U=UROOT .UTIP 
CALCs-DXLEOY ( J) /CHORD (U) 

CALC  =  CAlC-  (XCAP  (  I  )  *.5)  *r)CDY  (U)  /CHOHD(U) 
CALC=CALC*ChOhD ( J) 

CALC=CALC*DZnxu( I .U) 

IF(U.FU.U;3l.UR.U.F0.UB?.0R.J.EU.UTIP)  60  TO  170 
OZDYU (  I . J) =CALC* (ZU( I  * J)-ZU( I »U*1 ) )/ ( Y (U> -Y (U*l )  ) 
GO  TO  17? 

1  7n  uZOYi)  (  I  .U)  =CALC*  (ZU  (  I  .  J)  -ZU  ( I  .U-1 )  )  /  (  Y  ( U )  -  Y  ( U- 1 )  > 
IT?  CONTINUE 

00  18?  I=lLE.lTf 
DO  18?  U  =  JROUT » JT  IP 
CALC=-UXLtUY (J) /CHORD <U) 

CALCsCALC- ( XCAP ( I) *.5) *DCDY( J) /CHORD(U) 
CALC=CALC*CHONn (U) 

CALC=CALC*DZDXl ( I .U) 

IF  (U.EO.UHl  .0R.U.F(3.UB?.0».U.F,0.UT  IP)  GO  TO  180 
DZDYL  (  I  .U1 =CALC* (ZL ( I .U) -ZL (I.U*1))/(Y(U)-Y(U«1)) 
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GO  To  1B2 

lao  ozdyl(I»jj  =caLc*  (7L  )/(y(j)-y  (j-u  ) 

18?  CONTTNOt 

wWlTF  (6»<J49) 

00 

250  '^wire  (bt25i)  I.JHOOd) 
f^ETUPN 
END 


SOBHOUT  InE  ASPLINE(Nn»KO%Yr>,SO*(HO»XO»SO»YO»XP»YP»XPP»YPP»NFL) 

PAMAMETEWIZES  (XDtYD)  DATA  IN  TEWMS  OK  AHC-LENGTH  SO  ALONG  CURVE* 
THEN  SPLINE  FITS,  XP * YP * XPP ♦ YPP  AWE  DERIVATIVES  W.B.T,  SO. 

NKL.eO.l  FOP  UX/OSa  0  (INFINITE  OY/UX  AT  LEFT  ENDJ 

NKL.EO.a  FOR  l)Y/OSa*l  (POSITIVE  OY/DX  AT  LEFT  END) 

NFL. ED. .1  FOP  OY/OSs-1  (NEGATIVE  DY/DX  AT  LEFT  ENO) 


DIMENSION  XD(inO)  .YO(NO)  *XO(NO)  *YO(NO)  .Sn(NO)  .SO(NO) 
DIMENSION  XP(NO)  «YP(N0>  «XPP(NO) tYPPlNO) 


1 


EPSOsl.E-li) 

NlsNn-1 


SO(I)aO 

H1=0 

OXl=X0(2)-XD(l) 
OY1=YO(2)-Ylj(n 
ClsSOPT  <0X1**2*0YI*»2) 
SD(2)=C1 

IF (N0.EQ.2)PETUWN 
00  I  r=2*Nl 
OXl=XO (I)-XD (I-l) 
OYlaYD(n-YU(I-l) 
0X2aX0(l>l)-X0(I) 
aY2aY0(I*i)-YO(I» 
OXsXO(I*l)-XD(I-l) 
UYsYO(I-»l)-YD(I~l) 

CZaSORT (0X2**?^0Y2**2) 
C=S0«T(UX**2«’DY**2) 

A=  IDYI*0X-0X1*0Y) /? 
H=4*a/ (C*C1*C2) 

HAV=  (HWH)  /2 

0S=C1* ( 1 ♦ (C1/2«HAVI **2/6> 

SD(I)=SU(1-1)*0S 

CIsC2 


Hl=H 

CONTINUE 

DS»C1* ( 1+ (Cl/2*H) **2/6) 

SD(Nn) =S0 (NO-l ) ^DS 

CALL  5SPLINE ( ND *50 * XD »NO*SO » XO* XP » XPP » 1 ) 
CALL  SSPLlNE  (  '‘D*SO*YD.NO.SO*YO*YP*YPP*NFL> 
return 
END 


SUdPOUTiNt  SSPLlNE (NO* XU* YD ♦NO«X*Y»YP * YPP* NSW  ITCH) 

SPLINE  FITS  NO  DATA  POINTS  UD»YU)  .  ASSUMES  CUBIC  RIGHT  END. 
LOCAL  ARRAYS  0lY*n2Y«Q3Y  MUST  BE  APPROPRIATELY  DIMENSIONED. 

NSWITCH.EQ.I  FOR  CALL  FROM  ASPLINE  WITH  OX/OS*  0 

N6WirCH,E0.2  FOR  CALL  FROM  ASPLJNE  WITH  DY/DS«*1 

NSWITLH.ED.B  for  call  from  ASPLINE  with  OY/USs-1 

NSwITCH.EO.A  FOR  OIRECT  CALL  AND  CUBIC  LEFT  ENO 
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OiMfcNSlON  A|)(M))  .Y1)(^'U)  »X(HO)  (  'iO)  ,YP(NO)  ♦YPP(NO) 
DIMENSION  OlY  ( 1  iiU)  t))£’Y  nsr>)  .n3Y  (  r:)0) 

EPSI la-1 .t-lU 
EPSI?=-tPi>I  1 
,\iiMi=Na-i 
ux  =  xo  (<?)  -xo  (  1  ) 

If  <OX,£W.O. >  bO  ID 
0F=<Y0<2>-VU(  1)  )/{)X 
01 Y ( 1 ) a.b 

(SO  TO  (  l»c:t3*4)  ,NSWl  TCH 

naY  (  1  >  a3*  lOF-0) /ox 

00  10  6 

O^Y ( 1 ) =3*  <UF-1) /Ox 
GO  TO  6 

0?Y ( 1 ) =3* IDF  *1) /Ox 
GO  TO  b 
I?  =  3 

Ox1=XO(3)-XL)(£) 

DFla(YIJ(3)-YDl?)  )/l.)Xl 

C=  «DXl**2-0X**?) /OXl 

G=<0X^0X1 I ♦ (OX*?*OX 1 ) /nxi 

UlY (?) =C/B 

02Y (?) ab* (DFl-DF) /B 

DXaDX 1 

IJFaDF  1 

DO  7  IaI2tNIMl 
!)X1=XD  (  I  ♦  1  )  -XU  (  I  ) 

IF  (UXl  .EQ.O, )  GO  TO  36 
nFl=(YD(I*l)-tD(I)  )/(Ul 
Ba2* (OX*DXl ) 

FaG^*  (OF  1-UF  \ 

DENOMsB-DX*DlT  (  l-U 

n2T  ( I )  =  (F-l)X*D2Y  ( I-l  n  /OEMOM 

U1  Y  (1 ) SOXI/DENOM 

DXaOXl 

OFsOFl 

CONT INUE 

OX  I  a  XL)  ((YO- 1 1  -XO  (NU-2) 

CALCa  (0X*JX1  )  •  (0X^()X1*2)  /OXl 
CALCla (DXI**2-0X**2) /UXl 
OENUMs (CALC-0 lY (Nn-2) *CALC1 ) 

02 Y  (NO-l )  =  (F-02Y  (NO-2)  *CALCl  l/OtiYOM 


NIM23N0-2 
DO  H  I=ltNlM2 


K=NO-I- 1 

IF  (NSWlTCH.EO.^.&NO.iX.FQ.l)  GO  TO  0 
02 Y (K) =n2Y (K ) -01 Y (K ) *U2Y (K»l ) 

CONTINUE 

02Y(Nn)a( <OX*UXI)*02Y(NO-1)-DX*02Y(ND-2>  >/OXl 
IF  (NSv<  I  TCH.NE  .4)  GO  TO  9 
OXaxn (2) -XD ( 1 ) 

0XlaX0(3)-XD<2) 

02Y ( I ) a ( (UX^OXl) *02Y (2) -0X*D2Y (3) ) /DXl 
k;=no 

00  ll  IsUNIMl 
K  *  K  **  1 

oxi=xn(K*i)-xo(K) 

0FI=(YU<K+1)-Y0(K) >/UXl 

01V(K*l)aUFl'»OXl/6#(02Y(K)*2«O^Y(N4n  | 

0  3Y(K*l)a(D2Y(K4n  -l)2Y  (K)  )/DXl 
CONTINUE 

01 Y( n  =UFl-0Xl/6*<2*02y (I ) ♦U2V (2) ) 

03Y ( n  a03Y (2) 

IF (NSWITCH.EU, 1 )  GO  TO  16 

INTEPPOLAIING  Y 


00  IS  J=1«N0 
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r>on 


DO  12  I  =  l»Nl) 

DXsXn { I ) -X ( J) 

IF(UX,tiE.fcPSIl.ANn.DX.LE.tPSI2)  GO  TO  13 
IF(OX.GE,tPSI2)  r,o  TO  U 

12  continue 

60  TO  37 

13  y(J)=YU<I) 

YP ( J) =01Y  ( I ) 

YPP ( I ) =02Y ( I ) 

60  TO  15 

14  OX=X ( J) -XU ( I ) 

Y ( J) =YD ( 1 ) ♦DX* (DIY ( I ) ♦OX/2* (D2Y 1 1) ♦DX/3*0JY ( I ) ) ) 
YP  ( J)  =D1  Y  (I )  ♦UX*  (D2Y  <  I )  ♦DX/2*03Y  (in 
YPP ( J) SD2Y ( I ) ♦DX*03Y (1) 

15  CONTINUE 
60  TO  23 

INTEPPOLAT ING  X 

16  CONTINUE 

00  22  U=l»NO 
00  17  I=l*NO 

0Y  =  Y0(n-Y(Jl  ^ 

IF(0Y,6E.EPSIi.AN0.()Y,LE.EPSI2)  60  TO  IR 
IF (OY.GE.tPSId)  60  TO  IQ 

17  CONTINUE 
GO  TO  38 

18  Y<J)=YD(V 
X(J)=XD(I) 

YP ( J) sUl Y  < I ) 

YPP(J)=D2Y(I) 

GO  TO  22 

19  DXs-OY/OlY ( I )  ,  . 

20  YO=YO( I ) ♦UX* lUlY ( I ) ♦DX/2* (n2Y (I) ♦DX/3*D3Y ( I) ) ) 

IfToyTge^EPSI I .aN0«0Y.LE.EPSI2»  60  TO  21 
YOP»01Y  1 1 )  ♦DX*  (D2Y  ( I )  ♦0X/2*03Y  ID) 

OELXs-OY/YOP 
OXaOX^OELX 
60  TO  20 

21  X ( J) *X0( I ) ♦OX 

YP ( J) =01Y ( I ) ♦UX* (D2Y ( I ) ♦0X/2*D3Y (1) ) 

YPP ( J) s02Y ( 1 ) ♦0X*D3Y ( 1 ) 

22  CONTINUE 


CONTINUE 

WFTURN 

WHITE  (6.100) 

WRITE  (6.101)  XD(1).X0(2) 

STOP 

WRITE  (6.100) 

WHITE  (6.102)  I.XO(I)  .XOd*!) 

STOP 

WHITE  (6.100) 

4R1TE  (6.103)  J.X (J) .XO(Nn) 

STOP 

WRITE  (6.100) 

WHITE  (6.104)  J. Y (J) .YDINO) 

STOP 

FORMAT  (/5X.  IdHSUBHOUTINE  SSPLiNE/) 

format  (/bX.2lHERH0H  IN  INPUT  XU ( 1 ) =£ 1 2 , 4 .SX . 6HX0 (2) sE 1 2. 4/ ) 
format  (/5X.  16HERR0R  IN  INPUT  I  *  1 5 .5X . 6HX0 ( I ) =E 1 2.4 . bX ♦ 

[  8HX0( I*1)*E1?.4/)  ,  ^  ^ 

FORMAT  (/5X.2JHX(U)  IS  OUT  OF  HAN6E  J= 1 5. 5HX ( J) sE 1 2 .4 .5X » 

I  7mXD(NO)=E12.‘./) 

format  (/5X.23MY(J)  IS  OUT  OF  RANGE  J=I5.5X.5MY < J) *E12.4.5X. 
I  7HYD(Nn)=tl2.4/) 
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or;n  nnr^  o  Donoon 


r 


SUBMOUTINt  PROFL 


CALCULATES  WiNG/HODY  AW  GPIO  DATA  IN  THE  XI-ETA-ZETA  COMPUTATION 

domain. 


DIMENSION  ZETAL  <49. 19) »ZETAU (49. 19) 

COMMQN/CONST/  ALPHA « ZMACH .OPL IM. EPS I « wE . WO . BSP AN> A 1 2 ♦ CSA . I . flX2, 

►  SNA.rDETA.TOXI.TUZErA.OETA.OXr.DZETAflLE.IM.INOSe. 

'  ITAIL. ITE. JM.UROOT. JTIP.KM.KW. JBl. JB2» JH3.PLCB1. 

»  PLC02.PLCB3 

COMMON/COOHO/  CHOPOdS)  .OCOY  (19)  .DTOY  (19)  »0XLt0Y<l9)  .ETA  (19)  »E  (49) 
»  »FLE«EN.rTE.<;(l9)  »  GAMLE  .GAMN.GAMF  .6AMTE «  GAMkT  »H  (2S )  . 

‘  TAU(19) .X(49.19) ♦XCAP(49> , XLE ( 19) , Y ( 19) .ZETACZS) .6F» 

»  ZCAp(25>  .ZLEdR)  .0ZLEDYd9)  ♦ALPHATd9)  ,0A0Yd9)  .GWT 

COMMON/DA Tax/  A1 »A?.A3.A4.AS,ETA0»NtTA.NXl.NZtTA.XCAP0.Xl0»YCAP0. 

^  NF.XF  (ISO)  .ZFdBO)  .HFdBO)  .USECT  »  J0L*  JBT.NSU(5)  ♦ 

'  N5l,(5)  .YSIS)  »XLES(S)  .ZLES  (5)  .CS  (S)  »  ATS  (5)  *15(5). 

‘  FSTS)  .XSUdSO.5)  .XSLdSO.S)  .2SUd50.5>  .ZSLllSO.S) 

COMMON/SUHF  y  DELL (49.19) .DELS(49) ,0ELU(49»19) .DZDXL<49.19) . 

►  0ZDXU(49«19) .DZDYL(49.19) .0ZDYU{49. 19) »HrtBL(49.l9) « 

►  HWBU(49. 19)  .TZL  (19)  »IZUd9)  .0800(49)  .KBOOL  (49. 19)  ♦ 

‘  K800U(49.19) .ZL(49«l9) «ZU(49.19) 


304 

format 

1 

2 

306 

format 

309 

format 

1 

459 

format 

1 

461 

format 

601 

format 

1 

60? 

format 

dH-///lx2rtHSUwFAC£-ADJACENT  GHIO  POINTS/// 
-  -  *-  -  DOWN.  J  ACHOSS)// 


1A33HUPPEP  SUWF4CE  (I 
3X1H1 .3XlOHKRODUd  .J)  /) 

(I*..  (2016)  ) 

(///1X33HLOWFR  SURFACE 
lOHKHOOL  d . J) /) 
dH-///lX33HBTATI0NS  OF 
6HIZU(J) .9X6HiZL ( J) /» 

(lH-///lUoHNEGATiVF  UPPEK-SURFaCE  GhIO-POINT  OFFSET. 

24H  AT  PLANFOKM  station  I  =  I4«5H»  J  *  14/) 
dH-///lX40HNEGATrVE  LO»<tW-SU«F ACE  GRID-POINT  OFFSET. 
24H  at  PLANFORM  STATION  1  s  14, SM.  J  *  l4/) 


d  DOWN,  J  ACROSS) //3XIHI,3X, 
ZERO  STREAMWISe  SL0PE//3XIHU,9X, 


100 


120 

C 


200 

C 


270 


hpi=?.*atan  d . ) 
DO  100  1=1, IM 
DO  100  J=1,JM 
ZETAUd  »U)=0. 
ZETAL (I , J>  =0. 

DO  120  1=1, IM 

on  120  j=i,JM 

TC=TAU  <  J)  *CHORl) 
ZH=2u ( 1 , J) /TC 
ZeTAU< lfU)=ATAN 
Zh=ZL ( 1 ,U) /TC 
ZETAL ( 1 , J) =ATAN 


( J) 

( ZH/A5) 
(7H/AS) 


/HPl 
/HP  I 


AH=AS«HPI 
DO  200  1=1 
00  200  J=i 
ZH  =  ZETAU(  I 
TN1=T AN (HP 
HwBU (I . J) = 
Zh=ZFTAL  1  1 
TNlsTAN (HP 
HV.BL  (I  .  J)  = 


,1M 
,  JM 
,  J  1 
1»ZH) 

I ./ (AH* 
,  J  ) 
I*ZH) 

1 ,/ ( AH* 


(  1 
(  1 


♦TN1*TNI ) ) 
♦  fNl*TNn  ) 


00  270  1  =  1, IM 
DO  270  J=I.JM 
KHOUD  < I , J) =0. 
KHODL ( I .D) =0. 
DO  30  0  1  =  1.  I M 
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00  300  J=i»JM 
Zh=^ETAU(I»J) 

IF  (ZH.EQ.O.)  GO  TC 

00  «?80  ^=^W»KM 

If  (ZH.GE.Zfcla  (K )  ) 
KBOUIX  I » J) =K 
GO  TO  2ab 
CONTINUE 
2HsZETAL(  UJ) 

IF  (ZH.EO.O.)  GO  TO 
00  290  K=i»KW 
kkskw-k* 1 

IF  (ZH.LE.ZETA (^K) ) 

KROOL ( I « J) =KK 

GO  TO  295 

continue 

CONTINUE 

CONTINUE 

WWITE  (StJOA) 

OU  J05  I=I«1M 
wWITF  <b*J06)  r«(KF 
WHITE  (b»309) 

DO  310  1=1. IM 
WHITF  (6»30b)  I.(KE 


I.  (KftOnud.J) 

I  .  (KBOntd.J)  .Jsl 


DO  400  J=i»jM 
IZU(vJ)  =0. 

IZL  («J)  =0. 

NFX=NF-1 
UO  365  N=i.NFX 

IF  <PF(N.l) .LE.HF(N) )  GO  TO  370 

CONTINUE 

XF1=XF (N) 

00  375  l=lNOSt, ITAIL 
IF  /X(I.1> .GE.XFl)  GO  TO  300 

continue 

IFS=I 

JHOOTXssjRUOT-1 

DO  410  Jsl.jRUOTX 

IZU(J) sIFS 

IZL  ( J) slFb 

DO  450  JsJROOT.jTIP 

no  420  IslLE.jTE 

IF  (nZUXuU.J)  ,GT.O.)  GO  TO  420 

IZU ( J) =1 

GO  TO  425 

CONTINUE 

DO  430  IslLE.lTE 

IF  (OZOXLU.J)  .LT.O.)  go  to  430 

IZL  U)  =I 

GO  TO  450 

CONTINUE 

CONTINUE 

WRITE  lb. 4593 

00  460  Jsl.jM 

WRITE  (b.4613  J . I ZU ( J3 « I ZL ( J» 

DO  500  1=1. IM 
DO  500  Jsi.JM 
DELU(I.J)=0. 

DELLd.Dl'O. 

1)0  550  1  =  1. IM 
DO  550  J=1.JM 
KXsKRODU  d » J) 

IF  (KX.EU.O)  UO  TO  525 
OFLUll  tj>  =Z£  TA  (kX)  -ZETAUd  .  J) 

^X  =  KBODL  1 1 . J) 

IF  (KX.EO.O)  go  to  550 

DELL  (  I .  J)  =ZETALd  .  J) -Z£TA(KX) 

CONTINUE 

DO  bOO  1=1. IM 

UO  600  Jsl.JM 
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nnonno 


r 


IF  (f)ELU  (  1  ♦  J)  .Rf- .0.  )  GO  TO  b-JO 

wRITF  I,J 

STOP 

5'iO  IF  (HELL  (I  » JJ  .Ge  ,0.  )  GO  TO  GOO 
v^^^ITF  (6,002)  r,j 
STOP 

600  CONTINUE 
C 

DO  700  lalfIM 
700  DELS(I)=0»0 

00  710  1=IN0SE,ITAIL 
710  DELS(I)=GAMF 

DO  7?0  I=iLE,lTE 
720  DELS(I)=GAMWT 
IWX=IFS-1 
xFi=x ( iFs*  n 

ALCJ  =  ().S*  (  3.*YCAP0/ET  AO-hpi*a  4) 

ALC4=0.6* (mPI*A3-YCAP0/ETA0) / (tTAO*ET AU) 
yFS=PF ( 1 ) 

DO  722  N=2,NF 

IF  (YFS.Ll .MF (N) )  YFS=RF(N) 

722  CONTINUE 

DO  790  I=lNOSe,IFS 
J2=JflOD(I) 

J1=J2-1 
X1=X ( I , Jl) 

X2=X ( I , J2) 

Y1=Y ( Jl ) 

Y2  =  Y  < J2) 

SLOPE = ( YFS-Yl ) / ( Y2-Y1 ) 

XTEST  =  X 1 ♦SLOPE* (X2-X1 ) 

IF  (XTEST»GF,XF1)  go  TO  HOO 
XXXs( <XI-XF1)/XF1)**2 
YF laYFS*SUPT ( 1 .-XXX) 

XXXa ( (X2-XF 1 ) /XF 1 ) **2 
TF2sYFS*SUWT ( i ,-XXX) 

PAT lOa ( Y2-YI ) / (YF?-YFi) 

Yes(Yi-«Aria*TFi) /(I.-ratio) 

YCAPEsYE/bSPAN 

ETA1=ETA ( J2) 

730  ETAlSO=tTAl*eTAl 

FUNaFTAl* ( ALCJ^ALC4*FT A1SQ)-YCAPE 
FUNPRaALCJ^3.*ALC4*F.T  Also 
DEL=FUN/FUNPR 
erA2=ETAl-DEL 

IF  (DEL.LT . 1 ,E-10)  GO  TO  7S0 
ETA1=ETA2 
GO  TO  730 
7Sn  ETAE=ETA2 

GAM  =  ABS (E  TA ( J2) -FTAE) 

DELS ( I ) =GaM 
790  CONTINUE 
800  CONTINUE 
return 
END 


subroutine  coeff 

CALCULATES, FIXED  GEOMETRY-  AND  GHIO-RELATED  QUANTITIES  WHICH 
APPEAR  IN  The  transformed  POTENTIAL  EQUATION, 


COMMON/CAPS  /  CAPG(49,’9) ,CAP6T I ( 19,25) ,CAP6T2 <♦9, 19) ,CAPl ( 19) ♦ 

•  CAPH(*9  9)  ♦CAPHTU19,25)  ,CAPHT2(49,19)  fCAPiadQ)  f 

•  CAPGBTA  ’9) ,CAPIT (49,T9) 

COMMON/COEF  /  Al(49),».  49 )  ♦  A3  (49)  ,  AS  (49)  ,8 1  ( 1 9)  ,B2  <  19)  ,83  ( 19)  , 

•  B5(19) ,B> .19) ,C1 (25) ,C2(25) ,C3(25) ,C5(25) ,07(25) » 
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COMMON/CONST/ 


COMMON/COORO/ 


Dl»02*EltE2*F'l«F2 

ALPHAfZMACH.0PLl''<«ePSI»i«E»wG»BSPAN»AI2.CSA.RXl*RX2« 
SNA.TDETAf TOXl«tDZ£TA.DfcTA.OXl,D2ETA»iLt»lM»lNOSE# 
ITAIL* ITEt JMtJftOOT ♦ JT IP *KM . XK » JBl ♦ JB2t JB3»PLCbl * 
PLCB2»PLCH3 

CHOKD(19) .0C0Y(19»  *aTDY(19) .UXLEOY( 19) fETA ( 19) (49) 
tFLE»FN»FTE*G(19)  ♦6AMLE»GAM(^»6AMF»6AMTE»GAM*»T*M(?5)  * 


CHOKD(19) .0C0Y(19»  *aTDY(19) .UXLEOY( 19) fETA ( 19) (49) 
tFLE»FN»FTE*G(l9)  ♦6AMLE»GAM(^»6AMF»6AMTE»GAM*»T*M(25)  * 
TaU(19) ♦XI49«19).XCAP(49) tXLE(19) *Y ( 19) ,ZETA<25) ♦GF* 
ZCAP(25)  .ZLE<i9)  *UZ(.£DY(19)  tALPHAT(19)  »DA0Y(19)  .GwT 


PI=4,*ATAN ( 1 , ) 

ILX=ILE-1 
iTp=i  rE-»  I 
DU  20n 
CC=CHORD ( J) 

TT=rAU ( J) 
rc=TT*CC 

DXLY=0XLEUY (J) 

DZLYsOZLEUY ( J) 

OATy=OAoy ( J) *Pi/i8o. 

CIS=-OCOY (J)/CC 
CTS=-DTOy ( J) / fT 
CAPI ( J) =Ci5/CC 
CAPIR  ( J)  =  (C1S*CTS)  /TC. 

DO  loo  K=1,KM 
Z2=ZCAP (K) 

CAPUT  1 (J»K) =ZZ* (CIS^CTS) 

I  CAPHTl  ( J»K)  =2.*ZZ*  (CIS*CIS*CTS*Cri»^ClS*CTU) 

DO  120  1=1»IM 

CAPG( I . J) = (XCAP( I ) ♦0,5) *CIS-OXLY/CC 
I  CAPH(I,J)=2.*C1S*CAPG(I,J) 

ANGLE sALPH AT (U) *PI/lRO. 

TNAsTAN (ANGLE) 

SCAal./C05(ANGLEI 
00  140  I=1»ILX 
CAPGB(I.J)aO. 

CAPGT2(I»J)=-0ZLY/TC 
CAPHT2(I«J>s2,*CAPGT2(I ,J)*(CIS*CTS) 

So^UoMsIleIite 

XX*XCAP(I)>0.5 

CAPGB(ItJ>=TNA/TC 

CAP3T2( I » J)=-{U2LY-CC*XX*SCA*SCA*DATY*DXLY*TNA)/TC 
CAPHT2( I ♦ J) =-2,*SCA*SCA*0ATy*(UXLY-CC*XX*rNA*DATY) /rc 
1  ♦2.*CAPGT2(I,J)*(Clb^CTS) 

I  CAPIT ( If J) =SCA*SCA«DATY/TC*CAPG«( I tJ)* (CIS^CTS) 

00  lao  isiTPfiM 

CAPGR ( I f J) aO. 

CAPGT2 ( 1 f J) =- (nZLY-CC«SCA*SCA*UATyfDXLy*TNA)/TC 
CAPHT2 ( I f J) =2.*SCA*SCA*DATY* (OCUy  < J) ♦CC*TNA*DttT Y ) /TC 
I  ♦2,*CAPGT2(Ifs))«(CIS*CTS) 

I  CAPITdf J)aO. 

I  CONTINUE 

DXa2.*DXl*0Xl 
0E=2.*0ETA*0ETA 
0Z=2.*0ZErA*0ZeTA 
A1 (  i)aO. 

Al ( IM) *0. 

A2(i)a0, 

A2(lM)aO. 

A3( i)aO. 

A3( lM>sO. 

IMX=rM-i 
OU  210  I52»IMX 
Al(I)a(F(I*l)4F(I))/nx 
A3( I )a(E ( 1 ) ♦F ( I-l ) ) /nx 
A2<1)s-A1 (I)-A3(I) 


68 


AS (  1 ) =0. 

AS (^) =0  . 

00  ^20  I=3»IM 

?2n  AS( I ) = (F ( 1-1 ) ♦F ( T-2) ) /OX 
C 

Bi < 1 ) = (0 ( 1 ) *0 (2) ) /OF 
S3(l)=(b(l)*b(2) )/nF 
B2(1)=-S1(1)-H3(1) 

HI  ( JM) =0. 

H2 (OM) =0. 

JMX=JM-1 
DO  JIO  J=2«JMa 
HI (J)=<G(J^1) ♦G(J) )/nt 
d3  ( J)  =  (G(  J)  *0  ( J-n  )  /OE 
310  H2(J)=-rfl (J)-83(J) 
HS(1)=(G(2)>G(3) >/nE 
H5(2)=(G(1)*G(2) l/OE 
no  320  J=3»JM 

320  HS( J)  =  (6(J-1 ) ♦Gl J-2) ) /Oe 
JMX2=JM-2 
00  330  J=ltJMX2 
330  H7(J)=(G(J>1) ♦G(J^2) )/OE 
H7 ( JM-1 ) =U. 

H7(JM)=0. 

C 

Cl  (l)=u. 

Cl (KM) =0. 

C2( 1) =a. 

C2(KM)=0, 

C3(l)=0. 

C3  <KM) so. 

KMXsKM-I 
00  410  KS2»KMX 
Cl (K) s (M (K> 1 ) ♦H (K ) ) /Ol 
C3(K)s(H(K) ) )/0/ 
410  C2(K)s-Cl (K)-C3(K) 
CS(l)sO. 

C5(2)s0, 

00  420  KS3.KM 

420  C5 (K) 3 (H (K-1 ) {K-2) ) /OZ 
KMX2sKM-2 
00  430  Ksl ,KMX2 
430  C7 (K) 3 (H(K^1  )  ♦H (K*2)  ) /OZ 
C7 (KM-1 ) 3U  , 

C7 (KM) 30  . 

C 

01=1 ./ (4.*UXI«0ETA) 
02=-01 

El=l ./ (4.*0ETa*DZETA) 
E2=-F1 

Fl=l ,/ (4,*0XI*02ETA) 

F2=-F1 

return 

ENO 


c 

c 

c 

c 

c 


SUBROUTINE  INlT(MSTARf) 

initializes  solution  arrays. 


* 

» 


COMMON/CONST/ 


common/solvo/ 


alpha. ZMACHtORLI' 

SNA. TOE  TA. TOXI . Tl 


[M.EPSl  ..^E.WG.BSPAN.  AI2.CSA.RX1  .RX2* 
_  .  _  .  . .  ruZETA.DETA.OXl .U2ETA. JLE.IM.INOSE. 

I  TAIL . T  TE . JM, JROOT . JT IP.KM.Kw. JBl . JB2. JB3.PLCB1 . 
PLCH2.PLCB3 

OLl (49. IS) .0L2(49. 16) .DSl (49) .052(49) .OUl (49.16) ♦ 
0U2 (49. 16) . 0PM AX. OPMAXT .GAMMA (16) . IMAXI . ITERl . 

JMAXr  .  JMAXT.KMAXI  .K'1AXT.NSUP.P(49.  16.25)  .PLEl  (  16)  . 
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nor* 


*  PNtW(?5)  .PN0S£»Pi>YM(49»a5)  »PTU16.2S)  »PT2{16«25>  » 

*  PT3( lb»25) ♦PTtl (16) ,PTE2< 16) ♦PwRL (49*16) ♦PwBU(49*l6) 
C0MM0N/5UWF  /  DELL(49il6)  ♦ntLS(49)  ♦UELU  (4t'*  16)  *  OZOXL  ( 49  ♦  1 6 )  * 

*  0ZnxU(49*16)  *0ZI)YL(49,16)  ,OZUYU  (49*  16)  *HWBL(4V»16)  « 

*  HimR(J(49«  16)  *IZL(16>  *1ZU(16)  »JH00(49)  tKBOOL  (49*  16)  » 

*  KPOOU (49. 16) .ZL (49,16) .ZU(49*i6) 


100 

lOS 

iin 


120 

UO 

140 


145 

150 

300 


700 

750 

flOO 

900 


IF  (MSTahT.NE.O)  60  TO  105 
00  100  1=1 *1M 
00  100  0=1 »JM 
1)0  100  6=1. K4 
P( I , J.K) aO* 

1)0  110  0=1.  JM 
00  110  Ksl.KP 
Pri <J,K)aO, 

PT2  ( J.K)  a{j, 

PT3(J.K) aU. 

00  120  1=1. IM 
DO  120  Jsl.JH 
OLl (T,J)=0. 

0L2(I.J)=U. 

PWHL ( I » J) =0. 

OUl ( I . J) =0, 

0O2(  I  «  J)=l), 

PwBU ( I ♦ J) =0, 

DO  130  1=1. IM 
DSl (I)=0. 

OS2(I)=0. 

00  140  Jai.JM 
PLEl (J)=0. 

PTEl < J) =0. 

PTE2 ( J) =0. 

PNOSFaO. 

00  145  1  =  1, IM 
00  145  Ksl.KM 
PSYM< I .K) =0. 

IF  (MSfAPT.NE.O)  60  TO  300 
00  150  jal.JM 
6AMMA ( J) =0. 

HIa4.*Al  AIN  (  1  .  ) 
AN6LFxALPMA*P1 /1«0. 

CSA  =  C05(A'N6lE) 

SWAsS IN (ANGLE) 

PX1=1 ,/«t 
PX2=1 .-HXl 

rnxir2.*oxi 

T0ETA=2,«''ETA 
TDZETA=2,*0ZE( A 
AI2=1 ./ (2MACH«ZMACH) 

IF  (MSTaRI .EO.O)  GO  TO  900 
00  HOO  I=lNOSfc,IM 
JEX=JROO(I )-i 
DO  750  J=i»JEx 
KL  =  KROi)L  ( 1  ♦  J)  ♦! 

KU=KR00U( 1 .J) -1 
DO  TOO  KsKL.KU 
P( I ,J,K) aO, 

CONTINUE 

continue 

PETUWN 

END 


SUBROUTINE  SOLVE (ITER) 

8  EXECUTES  A  complete  RELAXATION  SWEEP  OVER  THE  COMPUTATION  DOMAIN, 
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nr>r>r>o 


C 


dimension  SLKWI <49.19) »SL2(19) .552 <19) »SU2(19) 
DIMENSION  SUKWl <49. 19) ♦SL3< 19) .SSJ(19) »SU3<19) 


COMMON/CAPS  /  CAP6<49.19) .CAPGTl (19.25) ,CAPGT2(49.19) ,CAPI (19) ♦ 
CAPH(49.19) .CAPHTl (19.25) .CAPHT2 (49. 1 9) .CAPIB(19) . 
CAP60(49.19) .CAPIT (49.19) 

ALPHA.ZMACH.0PLIM.EPSI.4E.WG.8SPAN.Ar2.CSA.RXi.RX2. 
SNA.TDETA.TDKI.TOZETA.DETA.OXI.OZETA^ILE.IM.INOSE. 


COMMON/CONST/' 


ITAIL.ITE.JM.JROOT.JTIP.KM.XW. JB1.JB2. J63.PLCei. 
PLCB2»PlC83 

COMMON/COORO/  CH0«0(19)  .DCOYdR)  .0TDY(19)  .0XLEDY(19)  .ETA(19)  .F(49) 


COMMON/SOLVO/ 


«FLe.rN.FrE.G(19)  .6AM(-£.6AMN.GAMr.GAMTe.6AM4T.H(2§)  . 
TAU(19).X{49.19) ♦XCAPT49) .XL? ( 19) . Y ( 19) ,2ET A ( 25) .GF . 
2CApt25) ♦ZLE<i9) .D2lE0Y(19) . ALPHAT ( 19) ,OAOY ( 19) ,6WT 
0Ll(49.p)  .0L2(49.19)  .QSl  (49)  .052(49)  .OUl  (49.19)  t 
002(49.19) ♦0PMAX.UPMaXT.6AMMA(19> .JMAXT .ITERT. 
JM4Xl.JMAXT.KMAXl.KMAXT.NSUP.P(49.l9.2$) .PLEI (19)  . 
PNEW(25) .PN0SE,PSYM(49.25) .PTl (19.25) ♦PT2(19,25) . 
PT3( 19.25) .PTEl (19) .PTE2(19) .PWBL(49.l9) .PW8UI49.19) 
COMMON/SUHF  /  OEl  L<A9.19) .DELS(A9) .0ELU(49.19).DZDXL(49.19). 

0Z0Xy(49.19) .DZDyL(49.19) .DZ0YU(49. 19) .HWBLlAR.lR) . 

1  I  It  U3  i  tlf  ^  1 • 


UPDATE  SUkFACE  CONDITION 

OPMAXsO.O 

NSUPsO 

ILXslLt-1 

IMX=IM-1 

INXaINOSE-1 

ITPalTE*! 

JMXsjM-l 
KMXsKM-  I 
OHs-DX  I 
OGs-GAMN 

CI»2.«0H*0H/ (OH *2. *06) 

C2»4,*0G/ (0H*2.*0G) 

C3= (nH-2,*0G) / (DH42,*UG) 

PW=P< iNUSt-l . 1 »KW) 

PWWsP ( I  NOSE -2. 1 .KW) 

PNUSF=C1*CH0RL’  ( 1 )  *C5A/FN*C2*PW*C4*PWW 
OGs-GAMLE 

CI*2.*UH*UH/ (uh.?,*OG) 

C2=A.*U6/  (L)H.2.*DG) 

C3s(l)H-2.*06)  /  (()H*?.*UG) 

CFsCSA/FLL 

00  20  JsUHOOT.JTIP 

PwaP ( ILE-1 » Jt^W) 

PW*»=P  (  lLt''2»  J.KW) 

20  PLEI (J)=Ci*CHUHn(J)*CF.C2*Pw.C3*PWW 

PLtl ( JHOOT-1 ) =2,*PLE1 (UPOOT)-PLtl (JkOOT.i) 

DO  30  HSUHFal.2 
DO  30  15EI=INUSF.IM 
JEssJHOD  ( ISET)  -1 
00  30  JSETsI.JE 

30  CALL  SURFbC ( ISET . JSFT .MSUhF) 

00  40  ISElalNOSF.IM 
JSETaJdOOdSET) 

40  CALL  SOPFSCdSET.JSFT.O) 

CALL  SYMBC(O.O.O) 

UHaOX  I 
i)6*GAM  T  L 

C1*2.*UH*UH/ (UG* (OH.OG) ) 
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ooo  ri  noo 


(OH-ob) ynG 
C3=  (OH-ufi)  /  ((;H*r>G) 

GO  ba  j=j^(ooTi jrrp 

PTE=(3.-0G/LiH)  *PWhL  ( I  TF-»  J)  -  1  1  .-OG/DH)  *PWHU  ( ITK-1  tJ) 

PE  =  P ( I  TZ*  1 » J»Kw) 

PEe=P  ( ITE^2* 

PTEl U) =Cl*Prt4C?*PK*C3*PFE 
50  PTE2(J)=3.*(PTEl(J)-PeUPee 

PTEl (JP00T-1)=2,*PT£1 (JPOOT)-PTEl (JROOT*!) 

COMPUTE  POTEMFIaLJ  REGION  UPSTREAM  OF  FUSELAGE  NOSE 

P ( INUSE  f i ♦««) =PN0SE 
IF  (AI2.GT,1,)  go  to  103 
00  102  J=1»JM 

no  102 

PT3(J»K)=P(ltJ»K) 
lOP  PT2 ( J.K) sP J,K) 

ISTART=3 
GO  TO  loS 

103  00  lOA  J=i.JM 
00  104  K=ifKM 

104  PT2 ( J»K) =P (  I  ♦ J«K) 

ISTAWT=2 

105  00  ISO  iSET=ISTART« INX 
00  120  JSET=1*JmX 

CALL  TRICUE <2«KMX. ISETfJSFT) 

CALL  invert (2«KMX) 

CALL  MAxH2tKMX4lSFT*JSETl 

00  US  KsPfKMX 

PTi (JSET,K)sP( ISFTf JSETfK) 

US  P  ( ISE  r»  JStT,K)  sPNFW  <K) 

PTl  (JSET«1)=P(ISET«JSET,U 
PTl ( JSET«KM) #P ( ISET.JSETtKM) 

120  continue  , 

00  125  KaUKM 
12s  PTl  (JM»K)aPUsET«JM«K» 

00  130  JsUjM 
00  130  KsifKM 
PT3< J»K>aPT2( J*n) 

130  PT2<J»IX)aPTl  (J«K) 

ISO  CONTINUE 

COMPUTE  potential:  REGION  OORNSTREAM  OF  FUSELAGE  NOSE 


DO  200  JslfJTlP 
200  P  ( ILF*  J»KR)  =PLEU  J) 

00  210  JsltJM 
SU2  < J) =PT2 ( J»^w) 

S03 (J) =PT3 ( J*KW) 
SL2<J>=PT2U*rtiy) 
SL3(J)sPTJ(J*<<;w) 

SS2 ( J) =PT2 ( J.NW) 

210  SS3<J)=PT3(J*i<w) 

PbAVEaP ( ILE* JROOT-1 tKW) 


00  son  1SET=IN0SF»IMX 

IS£TX  =  ISF,T-1 

IS£TP=IbFT*l 

JE2=jea0USETA>-l 

jeisJ8U0USET»-l 

JFOsJROOUbETP)-! 

IF  <  ISei  .EQ,  ILX)  OF.OaJEl 
IF  nSET.EQ.irP)  Jt2aj£l 


line  segments  below  WlNG/HOUr 

00  220  JsltJM 
PTP  < J*KW) abL2  < J> 

220  PT3 ( J*Kw) abL3 ( J) 

IF  (JE2.E'^.0)  GO  TO  230 
00  225  Jsi.jE2 
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c 

c 

c 


SUKWl  < ISETX, J) =P( ISFTX« J,KW*  1  ) 
K-^sKROOL  (  ISF.Ta  ,  j) 

P  (  ISFTX  ,  J»KH*  i )  =f)L  1  ( ISFTX,  J) 

7  IS  P  (  ISkTX»  =i)L2  (  ISF.TX,  J) 

230  continue 

1)0  23S  J=i*JEl 

SUKWl  (I5Er»j)=P(ISET»J»K*^+l) 
KB=KR0UL ( ISET , J) 
P(ISET»J*'^b*l)=nLl  (ISET.vl) 

23S  P( ISETf J«NH*2) =Dl2( ISET.J) 

DO  240  J=i»JEO 

5UK»*1  ( ISET P«  J)  =P  <  ISETP«  J»KW*  1  ) 

KH  =  KHO[)L  <  ISETP,  J) 

P  ( ISETPtOtKB*  i )  =I)L1  (ISETP.J) 

240  P(lSETP.J«KH*2)=nL2(lSETP.J) 

DO  25)0  JSL7  =  l,JEl 
KbsKROOL ( ISE  TtJSET) 

CALL  TPICUE (2»KH,1SET,JSET) 

CALL  invent (2*KH) 

CALL  MAXI <2.Kti, ISET* JSET) 

DO  245  Ks2«KR 

PTl (JSEI *K) =P( ISET, USE T«K) 

24S  P (ISET , JStTfK) =PNEw (K) 

PTl (jSET,i)=P(ISET,JS£T,l) 
K81=KB*  1 
00  247  K=Kril,NW 
247  PTl < JSET,N) =P I ISFT,JSET,K) 

CALL  SUPEdC ( ISET, JSET,  1 ) 

P ( ISET, JStl ,Kd* 1 ) =nLl ( ISET , JSET) 
P ( ISET, JStT ,Kb^2) =DL2( ISET, JSET) 
250  CONTINUE 

DO  255  J=1,JE1 
SL3 ( J) =SL2 ( J) 

255  SL2  ( J)  =PT1  < 

LINE  SEGMENTS  ABOVE  wiNG/HODY 


DU  260  J=i,JM 
PT2 ( J,KW) =SU2 ( J) 
PT3(J,Kw)=SUJ(J) 

IE  (JE2.E'-<,0)  GO  TO  7/0 
00  265  J=1,JE2 

P(ISETX,J,NW*1)=SUKW1 (1SETX,J) 
SLKwl  (  ISE I  X, J) =P ( ISETX, J,KW-1 ) 
K8=KB0DU ( ISETx , J) 
P(ISETX,J,KB-l)=nui (ISETX,J) 

P<  ISETX, J,Kb-2) =0U2( ISETX, J) 

CONTINUE 

DO  275  J=l,JEi 

P(lSET,J,^W4l)=SUKwl (ISET,J) 
SLKWI  (  ISET, J) sP ( ISET, J,KW-1 ) 
KB=KBOOU ( ISET , J) 
P(ISET,J,^8-1)=L}U1  (ISET,J) 
P{1SET,J,K8-2) =UU2(I5tT,J) 

DO  280  J=1,JE0 

P(ISETP,J,Kw*1)=5UKW1 (ISETP,J) 
SLK4  1  <  ISETP, J) =P ( ISETP, J,KW-1 ) 
KBsKflOOU ( ISETP, J) 

P  (  ISETP, J,NB-1 ) =001 (ISETP, J) 

P(  ISf  TP, J,Ke-2) =DU2( ISETP, J) 

DO  290  JSET=’ , JE 1 
K8=KP00U ( ISET, JSET) 

CALL  TRICuE (Kp, KMX, ISET, JSET ) 
CALL  INVERT (Kb, KMX) 

CALL  MAXI (KB, KMX, ISET, JSET) 

DO  285  K=KB,KMx 

PTl ( JSET ,K) =P( ISFT,JSET,K) 

P  (  ISET, JStT ,K) sPNEW (K) 

PTl ( JSE I ,KM) =P ( ISET, JSET, KM) 

KB1=KR-1 

DO  28 7  K=KW,KB1 
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or-.o 


?A1  pri  (  J'SF.T  ,K)  =P  (  ISFT»JSfcT»K) 

CALL  SUPFbC  (  IbF  r.JSfcT.*?) 

P(  1SFT»  JStT«Kd-l)  =(JU1  ( ISETf  JSET) 
P  ( ISf  T  t  J5LT  ,Kb-^)  =DU?  (ISETtJSED 
?9n  CO(v4TINUE 

DO  £95  J=l*JEi 
SU3( J) =SU£ ( J) 

£9*5  502  ( J)  =PT  1  ( JfKw) 

lE  (DE2.EO.0)  GO  TO  310 
DO  305  J=ltJE2 

305  P<  ISFTX,  JfKW-i )  =5LKI»1  ( ISETX«  J) 
31n  CONTINUE 

DO  315  J=l»JEi 

315  PdSET*  J»KW-1 )  sSLKWl  ( 1S6T,J) 

DO  320  J=i.JEU 

320  P(15FTP»J*Kw-1)=5lKw1 (ISFTPtJ) 


lines  outouapij  of  wing/body  side  edge 


325 
3325 

3330 

326 


3327 

327 

330 


33? 


335 


350 

370 
3  75 

400 

410 


P  < ILE« JWOOT-1 fKW) =P5AVE 
GO  To  3325 


IF  ( ISET.ED.  ILX) 

IF  ( ISET ,EQ.  I  rP) 

DO  325 
Pr2(J*Kw) =SS2(J) 

PT3( J.Kw) =SS3< J) 

GO  TO  326 

DU  3330  J=JROUT.JTIP 
PT2(J.K«i)=PTEl  (J) 

PT3 ( J«Kw) =PTE2 ( J) 
P(ITE.J»KW)=PTE1(J) 

P ( I TF-1 » J*KW) SPTE2 ( J) 

P ( I TE« JWOUT-1 ,KW) =PTF 1 ( JPOOT-1 ) 
IF  (JE2.EU.0)  GO  TO  327 
IF  ( ISET.LU. I  TP)  GO  TO  3327 
P  ( ISFTX,  jt2«Kv<)  =0S1  ( ISETX) 

IF  (JE2-l.LT,il  GO  TO  327 
P  ( ISFTX,  Jt2-1  ,K«;)  =nS?{  ISETX) 


CONTINUE 

PTl (JEl,K4)=nSl  (TSET) 

IF  (JEl.LT.l)  GO  TO  330 
PTl(JEl-l»KW)aUS2(ISET) 
CONTINUE 

IF  ( ISET.EU, ILX)  GO  TO  332 
P ( ISFTP, JtO,Krt) =0S1 ( ISETP) 

IF  (UEO-l.LT.n  GO  TO  332 
P  ( ISETP,U£0-1  ,Kw)  =r)S?(  ISETP) 
CONTINUE 
JSET=JBOf)(ISEf) 

CALL  SUKFbC ( ISET, JSET, 0) 
P(I5FT,JEI,KW)=US1  (ISET) 

IF  (JEl-l.LT.i)  GO  TO  335 
P( ISET, JEi-1 ,Ow) sDS? (ISET) 
JE=JE1*1 

(OO  370  JStT  =  JE,UBX 

CALL  TRICUE (2,KMX, ISET, JSET) 

CALL  INVEHT (2,KMX) 

CALL  MAXI (2, KMX, ISET, JSET) 

no  350  «=2»KMX 

PTl ( JSET ,K) =P ( ISET, JSET, K) 

P ( ISET, JSET, K) =PNEW (K) 

PTl (JSET, 1)=P( ISET, JSET, 1) 

PTl (JSET,KM)=P(ISET,JSET,KM) 

CONTINUE 

00  375  K=1,KM 

PTl (JM,K) =P ( ISET, JM,K) 

00  400  J=I,JM 
00  400  K=1,KM 
PT3(J,K)=PT2(J,K) 
PT2(«),K)=PT1  (J,K) 

00  ‘>10  J=1,JM 

SS3(J)=PT3(J,KW) 

5S2 ( J) =PT2 ( J,Kw) 
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r>oo  n  onono  nor-, 


DU  4?n  jsjEtJi'^ 

SL^  (  j)  =bS£;  (  J) 

4^0  SIJ2 ( j)  =bS£; ( j) 

500  CONTlNUt 

wRITf  (6,511)  I  fFR.OPMAX,  IMAX  I  »  Ji^AXl  »KMAX1  ,N5UR 
511  FOR'^AT  (IM  ,  I  iO,E  15,4, 315,  I  10) 

UPDATE  CIKCULAlION  AND  FAR-FIELD  CONDITION 

DU  oOO  JsJPOOT , JTIP 

5NFw=  ( 1?.-UG/1)P  )  *  (PWHU  (  I  TF  ♦  J)  -PW5L  (  ITE  ,  J)  )  “  (  1  .-OG/DH) 
1  (P^HU (ITt- 1 , J) -PwhL { ITE-1 , J)  ) 

600  GAMPA(J)=wG*6NF»i*(1,-wG)  *GAMMA  (  J) 

GAMMA ( JPOUT-1 ) =^,*GAMMA ( JPDOr ) -GAMMA ( JWOOT*l ) 

CALL  FAWHC 

^RlTf-  (6,011)  gamma  (  JWOOT)  ,  I  TEH  1  ,l)PMAXT  ,  JMAXT,KMAXT 
611  format  ( M*,55x,E15.4, I 10,E13.4,^I5) 

WE  turn 
END 


subroutine  FAHBC 

UPDATES  FAR-FIELD  BOUNDARY  CONDITION  IN  THE  TREFFTZ  PLANE. 


DIMENSION  SAVE(19) 
COMMON/CONST/  ALPH 


ALPHA,ZMACH,0PLIM,£PSI ,WE,W6,BSPAN,AI2,CSA»RX1,RX2» 
SNA,TOETA,TOXI,TUZ£rA,QETA,DXI,UZETA, ILEflMtlNOSEf 
ITAIL,  ITEt  JM,  JROOT  ,  JTIP,KM,K*(,  JBl,  JB2»  JB3,PLC8l» 
PLCB2«PLC83 


COMMON/SOLVO/ 


KM,K*(,  JBl,  JB2»  JB3,PL( 


OLl (49,19) ,0L2(49»19) ,051(49) »DS2(49) ,DU1 (49,19) ♦ 
0U2(49, 19) ,DPMAX»DPMAXT»6AMMA(19) flMAXl jITEBT, 

JMAXI , JMAXT.KMAXI »KMAXT»NSUP,P(49, 19,25) »PLE1 ( 19) , 
PNEW(25) ,PN0S£,PSYM(49,25) ,PT 1(19,25) ,PT2( 19,25) , 
PT3(  19,25) ,PTE1 (19) ,PTE2(19) ,PwBL (49, 19) ,PW3y (49, 19) 
DELL  (49,19)  , DELS  (49)  ,OELU  ( 4^^,  l9)  ’OZOXL  (49, 19)  , 
0ZDXU(49,19) ,OZDYL (49, 19) ♦0^0X0(49, 19) ,HWBL (49, 19) , 
HWBU (49, 19) , IZL ( 19) , IZU ( 19) ,JB0D(49) ,KBODL ( 49, 19) , 
K80DU(49,19) ,ZL (49, 19) ,ZU (49, 19) 


I tewt=o 

JMX=JM- 1 
KMX=KM-i 
JE=JBOD ( IM) -1 
100  DPMAXT=0.0 

00  110  J=l,Jt 

SAVE ( J) =P ( IM, J.KW.l  ) 

KB=KR0DL ( IM, J) 

P ( IM, J,KH* I ) =UL I ( IM, J) 

110  P ( IM, J,KB»2) =DL2 ( IM, J) 

DO  150  J=1,JE 
KH=KBOOL ( IM, J) 

JSET=J 

call  TWIT  (2,Kti,0,jSKT) 

CALL  INVEWT(2,KP) 

DO  125  ^=2,^ri 
OFLP=PNEw(K)-P(IM,j,K) 

()P  =  AP5  (DELP) 

IF  (DP.LT .DPMuXT )  GO  10  125 

dpmaxt=dp 

JMAXT=J 

kmaxt=k 

PT  1  (  J,K  )  =p  (  IM,  J,k:  ) 
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rioorvi 


•f 


I3n 

150 


2ln 


225 


230 

250 


310 


P ( IM, J,K ) =HNE«  (K ) 

MT 1 ( J» 1 ) =H I IMt J, 1 ) 

KBl=KB>i 
no  130  ^s^Hl»^W 
P71 (J^K)=P(IM,J^K) 

CONT INJE 

00  210  Jrl,jE 
P  { IM,  J,Kw*l)  =5AVF,  ( J) 

S&Vt ( J) =P ( 1M*0»KW-1 ) 

KB=KP0D0(iM, J) 

P ( IM, J»KB-1 ) =UU1 ( IM. J) 

P  (  IM,  j»KB“2)  =0U2'UP«  J» 
on  250  J=l* JE 
KBsKROOU ( IM, J> 

JStT  =  J 

CALL  T«IT <KB»^MX,0»JStT) 

CALL  INVENT (^B,KMX) 

DO  225  K=^B,KMX 
OELP=PNEw(K)-P(IM,j,K) 

OPsARS (HELP) 

IE  (OP,LT,UFMAxT I  GO  ro  225 

opmaxt=op 

JMAXT=J 
KMAXT  =K 

PTl ( J,KI =P( IM, J,K) 

PnM,J,K)=PNE«  (K) 

PT  1  < J,KM) =P ( IM, j,km) 

KB1=KB-1 
00  230  K=KW,Kr!l 
PTl  ( J,K) =P ( IM, J,K) 

CONTINUE 

00  310  J=i»JE 
PaM,J,KW-l)=SAVE(J) 
P(IM,jE,KW)=DS1 (IM) 

IE  (JE.NE.l)  P(IM,vJE-UKW)sOS?(IM) 
PTl  (je.KW)sOSUIM) 

IE  (JE.NE.l)  PT1(JE-1*KW)=0S2(IM) 

00  350  JaOHOOT ,JMX 

JSETsJ 

CALL  TBIT  (2,KMX,0,OSC T) 

CALL  INVENT |2»KMX) 

00  325  k=2,kMA 
OELPsPNEP (K) -P ( IM,o,K) 

OP=ABS (DELP) 

IE  (OP.LT.OPMAXT)  go  to  325 

0PMAXT=DP 

JMAXTaJ 

HCMAXTaK 

PTl ( J,K) =H ( IM, J,K) 

P ( 1»,J*K) (K) 

PTl  (J,l)=P(IM,j,l) 

PTl ( J,KM) =P ( IM, J,KM) 

CONTINUE 

iteht=iteht*i 

IF  (  ITEPT.Ot  ,<iOO)  PETuPN 

IE  (HPMAXI.bT.UPLIM)  GO  TO  100 

PETUPN 

END 


SUSPOUTINE  SURFflC ( I « J»M) 

CALCULATES  INTERIOR  IMaGE-POINT  POTENTIAL  VALUES  REOUIREO  TO 

Satisfy  the  flow-tangency  boundary  condition  on  lower  (m,eo.i> 

AND  upper  <M,EQ,2)  SURFACES  OF  THE  WIN6/B00Y, 
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c 


COMMON/CAHs  /  CAPG(49,19i . C APGT 1 ( I 9 , 2b) .CARGlZ (49, 19) ,CAPI (19) . 

COMMON/CONST/  §fc;H?5i«iv6fc(3ilTSi5in:§;!?£J^TJ!f;HtA:W6?ff’ 

‘  ITAIL* ITE» JM, JROOT, JTIP,KM,KW, J81 » J82» JB3,PLCB1 , 

CH0HDjl9) ,DCDY<19) .UTOY (19) .UXL£DY{19) ,ETA ( 19) »F (49) 
jFLE jFN»FTE«6( 19) ♦0AMLE»6AMN*6AMF»GAMTE»GAMiT,H(2b) ♦ 

«XlE(19)  ,Y(19)  ♦ZETA(25)  tGF, 
2CAP (25) ,ZLE (19) ♦UZlEUY (19) # ALPHA! (19) ,OADY (19) ,691 
OLl (49, 19) fOL2(49»  19)  »OSl  (49) ,0S2(49) »OUl (49, 19)  » 

OU2 (49, 19) ,OPMAXjDPMaxT,GAMMA ( 1 91 , I  MAXI , I TEPf » 

JMAXI , JMAXT,KMAXl »KMAXT ♦NSUPtP (49, 1 9,25) ,PLE1 ( 19) , 

85bt«|5!n!';S!!§(f?U?!iV(S?a«U8Jf5V.'J28i?U.,„, 

MWBU (49, 19),12L(19),IZU(19) , JB00(49) ,K800L (49, 19), 
KB0DU(49,19) ,2L ( 49, 1 9) , ZU (49, 1 9) 


COMMON/COOHD/ 


COMMON/SOLVO/ 


COMMON/SUHF  / 


()k  =  (JE7A 

rc=T&u  (  j)  *CHOr(o  ( j) 

IF  (M.eQ.2)  GO  TO  100 
KP  =  -l 
UL=“0Z£TA 
DEL=-0£LL ( I , J) 

IZ=IZL ( J) 
bLOPX  =  OZOXL (I , J) 

SLOPY=UZDYL (I » J) 

KS=^P0DL ( I » J) 

H«»H  =  Hv^8L  (  i ,  J) 

CAPGTH* (OL-DEL) * (CAPGT 1 ( J,Kb) -CAPGTl ( J,KS*1 ) ) /DL*CAP6T1 ( J,KS* 1 ) 
C4PGTe=CAPGT8*CAPGT2 ( I , J) 

GO  TO  120 
100  KP=,1 

OLsDZETA 
()EL  =  DELU  ( 1  ,  J) 

IZsIZO(J) 

SLUPX=OZOXU (I , J) 

SLOPysOZDYU ( I » J) 

KS  =  KBODiJ(  I ,  J) 

H,rB  =  H»K80(  I  ,J) 

CAPGTR= (DL-DEL) o (( flPGT 1 ( J , KS ) -C aPGT 1 ( J,KS-1 ) ) /DL*CAPGT1 { J,KS-1) 
CAP6TO  =  CAHGTt3*CAPGT2  (  1,0) 

120  IF  (I.EU,iM)  GO  TO  210 

IF  (J.Gt.OPQOr)  GO  TO  140 
aH=-nx  I 
IP=-1  , 

bO  TO  200 

140  IF  d.GE.IZ)  GO  TO  160 
OM  =  -l)X  I 
IP  =  -1 
GO  TO  200 
160  DH=,nXl 
IP=M 

200  CaPK= (bLOPX/CMUWO (0) ♦CAPG ( I ,0) *bLOPY) ( I ) 

CAPL=SLOpr*G( J) /R5PAN 

CAP’^sHtVb*  (CaPGR  (  I  ,0>  *bLOPX,CAPG(H*SLOPY-l  ./TC) 
CAPP=CSA«bLOPX-SNA 
GO  To  21b 
210  CAPK=0.0 

CAPL=SLUPY*G(0) /hspan 
CAPMrHWri* (CaPGTH*SLOPY-1 ./TO 
CAPPs-SNA 
GO  TO  220 
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21S  Al=-1.5/i)n 
A2  =  ^.f)/DH 
A3=-0.S/0M 
i>0  Tf)  £?^b 
220  Al=0.0 
A2  =  0.0 
A3=0.0 

2ZS  t3l=-1.5/0^ 
tl2  =  2.i)/UK 
B3=-0.b/DK 

Cl=-  (nL*2.*L)EL)  /  (DtL«  (0L*l)EL>  ) 
C2=(r)L-‘DEL)/(L)FL*DL) 

C3  =  -nEL/  (l)L*  ( JL^OFL  )  ) 

01  =  0.S*  (DL*nEL»  *  (2,*DL^nF.L)  /  (DL*OL) 
i)?=-nEL«  (2.«0L*nF,L)  /  (nL*r)L) 
03=u.b*UEL<»  (UL*OEL»  /  <L)L*OL) 
£i=2.«dl*ul/  (ut.L<*  (ni^uen  > 

E2  =  -2.'^(DL-DEI.)  /OFL 
E3= (OL-UEL) / (UL^OFL) 


Jb  =  J*  1 
Jfi  =  J*2 
KS 1 =KS*^P 
Kb2=KSl ♦KP 
P1=H ( I » J.^S) 

P2  =  P( I  * J«0S1 ) 

P5=D1*P ( I > J5»^S) ♦D2*P (1 » Jb«KSl > ♦Oi*? (I . Jb»KS2) 
P6  =  01«P  ( I » Jb^^SI  ♦OP<»P  ( I  »  J6.KS1 )  ■*03*P  ( I .  J6.KS2) 
IF  ( J.EQ.JHl-l  ,0P.  J.EiJ.JRl »  GO  TO  181 
IF  (J.EO. JB2-1.0R, J.E0,JB2)  60  fO  1182 
IF  ( J.EO. 0B3-1  .OP. J.EO, 0*33)  60  TO  182 
GO  TO  118^5 

181  JKsJRl 
PLC=PLCBl 
GO  TO  183 

lie?  jKsjH? 

PLCaPLCB2 
GO  TO  183 

182  JK=J83 
PLC=PLtb3 

183  CONTINUE 

PA=(2.*P  (1 1 JK*  1  ,KS)  ♦  <  1  .-PLC)  *(E»n  ♦  JK-1  tKS) 

1  -2.*P( I^J^f^5) ) )/(l.*PLC) 

PB=  (2.*P  ( 1 » JK^l  .Ksl  )  ♦  ( 1  ,-PLC)  *TPT  1 » JK-1  »^Sn 
1  -2.*P(IiJK»KS1) ) )/(l.^PLC) 

PC  =  i  *  J'<-*  1  fKS2)  ♦  ( 1  ,-PLC)  *  IP  ( I  ♦  J^-l»^b2) 

1  -2.*P(I»JK,K82) ) )/(l.*PLC) 

IF  (J.EQ.Ji^)  GO  TO  184 
Pft=ai«Pd*02*PB*03*PC 
GO  TO  118b 

184  P5=01*PA*U2*P8^03*PC 

PAA  =  P  (  I  f  J^•'l  »KS)  ♦3.*(PA-P(I  ♦  JK»KS>  ) 

PeB  =  P  ( I  f  J^-l  tKSl )  ♦S.*  (PB-P(i  f  OKtl^Si)  ) 

PCC  =  P ( I » J^-l »^S?) ♦B.* (PC-P( I ♦ JK»^S2)  ) 

P6=0l *PAA^D2*PBB*D3*PCC 
lies  IF  (I.Nt.lM)  GO  TO  230 
P3=U,0 
P4s0,0 
60  TO  240 
230  I3=I»IP 
I4=I3>1P 

P 3  =  01  •PC  li.jfKS)  ♦n2<»P<I3»J»KSl>  ♦r)3*P(I3»J»KS2) 
R4  =  01*P  (  I4,J.KS)  ♦D2*P  ( I4»J»KS1  >  ♦l)3*P  ( I4»  J»KS2) 
240  OENOM=CAPN*A1*CAPL*H1 ♦CAPM*C1 
T 1=-CAPK* ( A2*P3*A3*P4) 

T2=-CAPL* (rt?«Pb*B3*P6) 

T3=-CAPM* (C2*Pl^C3*PP) -CAPP 
PBUPFs ( T l*T2*r3) /DENOM 
P01=F l*PSUWF*E2*PI*F  3*P2 


P01=F l*PSUWF*E2*Pl*F  3< 
PD2=3.*PD1-3.*P1 ♦PP 
IF  (M.EU.P)  GO  TO  250 
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OLI < I ♦ J) =H01 
UL^< 

PWPL  (  I  *  J)  =PSlJr(F 
WETURN 

PSO  OUl  (  I  ♦  J)  =PI)1 
OU^ ( I . J) =PU^ 
RWBU ( I » J) =PSUHE 
RETURN 


PLANEORM-bIOE-El)fjE  1  N'Aut -Pdl  NT  POfENTlALS 

ENTRY  bURESC 
IMa=IM-1 
UK  =  L)E  TA 
0EL=0e:LS  ( i ) 

IF  (I.GT.ITE)  GU  TO  100 
DH  =  -f)Xl 
IP  =  -1 
GU  TO  JIO 
300  0H=*nXl 
IP=*1 

310  FP  =  F  (  I  ) 

IF  (  I  .LT. ILE.UR. I .OT.  ITE)  GP  =  GF 

IF  (1 .GE.ILE.AND.I.LE.I  TF)  GPsUwT 

GSET  =  (  (UK-OE'L)  *CAPG<  I » J)  ♦DFL*CAPG(  I  ,  J-1 )  )  /DK 

GTA=CAPGT1 < J.Kw) ♦CAPGT2 ( I . J) 

GTB=CAPGT1 ( J-1 «KR) ♦CAPGT2 ( I ♦ J-1 ) 

GTS£T=( (U^-OFL) oGTA*DEL*GTB) /D^ 

Al=l .5/JH 

A2=-2.0/OH 

A3=0.5/DH 

Rl=- (UK42.*D£L) / (DEL* (Ok^OEL) ) 

02= (DK^UEL) / (uEL*0K) 

03=-OeL/ (UK« (OK^OEL) ) 

£1=2,«0K*JK/  (JEL*(r)K*DEt)  ) 
E2=-?.*<DK-DEL)/DEL 
£3=<DK-UEL) / (UK ♦DEL) 


12  =  1  IMP 

IF  (I.EQ.IM)  llrl 
IF  (UEU.lMx.OH.I.EO.IM)  I2=Il 
Ji=J^  I 

P1=P  (  I  *  J*N**) 

P2  =  P  ( I  *  J1  -K**) 

P3=( <D^*OtL)*P{IltJ»KW)-OEL«P(Il»Jl»KW) )/UK 
P4= ( (OK^OEL) *P ( !?♦ J»Kw) -OEL*P ( 12» Jl tKw) ) /UK 
Pb= ( (OK ♦DEL) *0 ( 1 ♦ J.Kw-1 ) -0EL*P ( I » J1 fKw-l ) ) /DK 
P6=( (0K^uEL)*P(I»J*KW^l)-DEL*P(I»JliKW4l) ) /OK 
GG= ( (DK^OtL) »GAMMA ( J) -nEL*6AMMA ( J1 ) ) /UK 
IF  (I.GT.ITE)  P5=P0*GG 

SCl=FP*GStT 
5C2=GP/HSPAN 
SC3=H(KW) *GTSET 
IF  (I.EU.iM)  '3C1=0. 

0EN0M  =  SC1<»A1>SC2*H2 
T1=-SC1* { A2*P3*A3*P4) 

T2=-SC2* (b2*Pi4B3*P2) 

T3  =  -SC3* (Ph-Pb) / (2.*0/ET  A) 

PSURF= ( T 1*T2* I  3) /DFNOM 
USl  (  I  >  =E1«PSURF^E2*PUE3*P? 

0S2(n=3,*0Sl  (n-3.*Pl*P2 
RETURN 


WING/HOUY  bYMMETRY-PLANE  POTENTIALS 
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41f) 


420 


470 


480 


490 

500 


ENT^iY  SYMWC 
KWX=KW-1 
KngPsKW*  I 
Ki4X  =  KM-  I 
00  500  II=2«1^ 
r,SET=CAPGi  1 1 «  1) 

FCl=GSfcT«F  (ID  /IJXI 
FC2  =  G(l)/(BSP4N*t0ETA) 

KLsKHOOL ( 11 « 1 > 

KUsKHOOU ( 11 ♦ D 
IF  (II.EO.IM)  GO  TO  470 
00  410  K  =  2«Ki/JA 

IF  (KL.NF.O.ANO.K.GT.KL)  go  to  410 
CAPGT=CAPGTl ( l,K) ♦CAP6T2(1I»1) 

FC3=CAPGT*n (K) /D7FTA 

IF  (GSET.LE.O.)  PSYM (I  I ♦K )  =  <FC 1 ♦PS YM ( 1 1 -1 »^ ) -F C2*P < 1 1 » 1 *K ) 

1  ♦FCB^PSYM (I  I tK-1 ) ) / <FC1-FC2^FC3) 

IF  (GSET.oT.O. )  P5YM<II»K)  =  (-FC1*PSYM(ID1»K)-FC2^P(I1»J»'^) 

1  ♦F  C3*P5YM ( 1 1 .K-1 ) )/ (-FC1-FC2*FC3) 

CONTINUE 

DO  4?0  K=NWP«KMX 

tXKsKMX-K^Kwp 

IF  (KU.NE.O, AND. KK.lt, KU)  GO  TO  420 
CAPGTsCAPGTl  (DkK)  ♦CAPGT2  ( 1 1 , 1 ) 

FC3  =  CAPGT*H(KK) /nzF.TA  .  . 

IF  (GSET.LE.O. »  PSYM  (  H  »KK )  s  (FC  1  *HSYM  (  H -1  tKK  ) -Fr.2*P  (  1 1  ,  1  tKK ) 

1  -FC3*PSYM(  1 1  .KK^D  )/(FCl-FC2“FC3) 

IF  (GSET.GT.O.)  PSYM ( 1 1 ,KK >  =  ( -FC 1 •PbYM ( 1 1 ♦ 1 » KK ) -FC2*P ( 1 1 ♦ 1 1 kK ) 

1  -FC3^PSYM(II,kk»1) ) / (-FC1-FC2-FC3) 

CONTINUE 

IF  ( 1 1 .6E. INOSE)  GO  TO  500 
CAPGTsCAPGTI  (  l,KW)  ♦CAP6T2(n,l) 

FC3=0.5*CAP6T*H{KW) /PZETA  _  ,  ,, 

PSYM(II,Kw)s<FCl»PSYM(Il-DKW)-FC2*P(II»l»KW)-FC3*lPSYM(Il,KW*D 
1  -PSYM(IIfKW-D  M/(FC1-FC2> 

GO  TO  500 
DO  400  Ka2,KL 

CAPGTaCAPGTl ( 1 ♦ K ) AP6T2 ( 1 1 « I > 

FC3=C4PGT*H(K) /DZETA 

PSYM( IM,K)s(-FC2*P( IM« l,K) ♦F C3*P5YM ( IM,K-1 ) ) / (-FC2*FC3) 

00  490  K=KU»KMx 
KKsKMX-KoKU 

CAPGT=CAPGT1 (i,KK) ♦CAPGT2(II,1) 

FC3=CAPGT*H(Ka)/DZE7A 

PSYM ( IM,KK) = (-FC?*P ( IM, 1 ,KK) -FC3*PSYM ( IM^KK^D ) / (-FC2-FC3) 
CONTINUE 

return 

END 


SUeROUTjNE  TRIC0E(K1,K2,I,J) 

CALCULATES  FINUE-OIFFERENCE  COEFFICIENTS  AT  GRID  POINTS  ALONG 
The  segment  ki  to  k2  of  line  (I»j). 


JgHH 


MMON/CA 


COMMON/COEF  / 


CAPH(49»t9)  fCAPHTl  (l9*2S)  »CAPHT2(49t  19)  fCAPiBh9)  ♦ 
CAPGB<49*19) fCAPiT I49«T9) 

Al  (49)  f  A2(49)  *A3U9)  *A5(49)  »Bl  (19)  ,82(19)  ,83(19)  , 
85(19) ,97(19) ,C1 (25) *02(25) ,C3 (25) ,C5 (25) ,C7 (25) ♦ 
D1,D2,E1,E2,F1,F2 
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10 


20 


30 


COMMON/CONST/  ALPHAt/rMACH,DPLlM»E.PSI  tWEtW^f  BSPAN 
SNA.TDETA.TDXI .TDZElAtDETAtDXlfOZE 
ITAIL»ITE«  JM*  JHOOT*  JTIP»KM.f<w«  JBl  ♦ 
PLCH2.PLC83 

COMMON/COUHD/  CH0PD<19) .OCDY(iy) *0TDY(19) tDXLEDY 
»FLE»FN»FTE«6( 19) »GAMLE»GAMN»GAMF» 
TAU(19) «X(49*19) ♦XCAP(49) «XLE(19) « 
2CAP (25) tZLE (19) tOZLEDY (19) ♦ ALPHA! 

coHHON/soLvo/ 

JMAXI ♦ JMAXT.KMAXl*KMAXT»NSUPtP(49» 


«Al2fCSA*PXl*PX2» 

TA«lLEtIM«INOSE* 

JB2*JB3*PLCB1» 


PNEW(25)  .PN0SE»PSYH(49.25)  tPTJL  (19f 
PT3 (19.25) .PTE  1 (19) »PTE2(19) .PW0L( 


(19) .ETA(19) .F(49) 
GAMTE.GAMwT,M(2b) . 
Y(19) »ZETA(25) .GF, 
(19) ,DA0Y(19) .6WT 

tSlJjVHIIfi’'* 

19*25) «PLE1 (19) « 
25) .PT2(19.25) * 
49*19) *PWBU(49*19) 


CHECK  P.D.E.  type 


FPsF (I) 

GP»G( J) 
CJ«CHORO(J) 
TC*TAU(J)*CJ 
C0G»CAPG(i*J) 
CU1»FP/CJ 
CU2-*CAPGB  ( I  .  J) 
-V|»C9G*FP 


CV2=GP/0SPAN 


^A^Iafp/ KJ*CJ) 


^-»FP*CfaG«CaG 
CAP3=-2.*FP*CdG/CJ 
C8P1*GP/ (BSPAN*BSPAN) 
CCP1*CU2«CU2 
CCP2*1 ./ ( TC*TC) 
CCP3*-2.*CU2 
CCP4a-2./TC  , 
CCP5a-2.*CU2/rC 
CDPla2.»FP*6P»CBG/BSPAN 
CDP2a-2.*FP*QP/ (CJ«BSPAN) 
CEPla2.*6P/BSPAN 
CEP2a-2,*GP*C02/BSPAN 
CEP3a-2,*GP/ ( rc*RSPAN) 
CFPla2.*FP*CU2/CJ 
CFP2a2i*FP*C8G 
CFP3a-2,*FP/CJ 
CFP4a-2,»FP*CU2*CBG 
CFP5a-.2,*FP*CdG/TC 
CFP6a-2,*FP/ (TC*CJ) 

■  '  - ( J) *CJ 


'jPla2,*CAP 

'JP2=-CAPH( 

JP3a2,*rAP 


AP 


^JP4a2;* 
CJP5a-TC 
CJP6a-CU2*TC 


.J)*Cv 

T jif j; 


0  ( J)  * 


*TC 

C 


JKaJM 
JKPaJK* 1 

IF  ( J.EU.JBl .OR.J.EQ.JBl*! ) 

IF  (J.EQ.JB2.0R.J,EQ.JB2*1) 

IF  (J.EO.JB3.OR.J,EO.J03*i) 

GO  TO  90 

JKaJBl 

PLCaPLCBl 

GO  TO  00 

JKaJBZ 

PLCaPLC02 

GO  TO  00 

JKaJ83 

PLC»PLCB3 


GO  TO 
60  TO 
GO  10 


10 

20 

30 


81 


noon 


80  JKH=JK*1 
90  CONTINUE 

00  700  K=K1,K2 
HP*H»K) 

GTP=CAPST1 (J»K) ♦CaPGT2( I*J) 

HTP=CAPHTI(J»K) ♦CAPHT2(I,J) 

Pi  =  P( 

P2*P( 

P3»PT2{J,5) 

p5spn»j»K-n 

P0»P|l»J»K*li 

IF  (J.EQ.l)  PU»PsyM(I*K> 

IF  (J.NE.i)  PilaPTl <J-1*K) 

P16»P(I»J>1»K) 

IF  (J.EQ.JK)  PI6=(?.*PI6^(1.-PLC)*<P11-2,'»P2)  >/(l.*PLC) 

IF  (J.EQ.JKP)  P11=<PLC*P16-2.«<PLC*P2-Pll) )/{2.-PLCJ 
UP=CSA*CUl*(Pi-P3) /T0XI*CU2*HP*tP8*P5)/TDZ£TA 

VPsCVl* (P1-P3) /T0XI*CV2*(P16-PIH /TOETA^GIP^HP* (P8-P5) /TDZETA 
WP*SNA*CWl*HP* (Pfl-PS) /TOZETA 

{lll^lli^wfiNS?KlNEi^§>U  GO  TO  105 
UP=UP-CU2*hP*GAMMA ( J) /TDZETA 
VP*y/P-GTP*HP»GAMMA  ( J)  /TOZETA 
WP*WP-CW1*HP*GAMMA ( J) /TOZETA 
105  UPSQ*UP«UP 
VPS(lsVP*VP 
WPSU*WP*WP 
(}PSQ«UPSQ*VPSO*WPSO 
TQsAr2«0,20-l .20*QPSQ 

apso»to*gipsq 
IF  (TQ.LT.O.)  go  to  200 

DEFINE  POINTS  OF  ELLIPTIC  (SUBSONIC)  COMPUTATION  MOLECULE  AND 
COMPUTE  TKI-OIAGONAL  COEFFICIENTS 

P3NaP(I-l»JfK) 

IF*n*!Eilt*TLANO.Kl.EO.KW*l)  P4*P44GAMMA  ( J) 

P6Nap(I-l»U«K-l) 

P7*P(1*1»J,K*1) 

P9N=P<  I-l  ♦J.K'H) 

P15=P<1*1*J*1«K) 

Pl7NaP(I-l»J*i,K) 

P18ap(I,J*i,K-n 
PI9aP(I,J*I,K*l) 

PllNaP(I.J-I,K) 

P12NaP(i-l,U-i,K) 

P13Nap(I,J-l ,K-1) 

P14NaPn.J-l»K»l) 

IF  (J.EQ.UK)  GO  To  112 

GO  18*125'^'^^* 

110  PlO»PSyM{l4l fK) 

PllNaPSYMd.K) 

PlgNaPSYM  I-1*K) 

Pl3NaPSYM(l,K-l5 
Pl4NBpSyM<IfK^t) 

GO  TO 
P15*(2 

PTNaJ. 

P  9a(2,*P19* (I.-PLC)*<P14N-2,*P8) )/ 

60  TO  120 

P10a<PLC*P15-2.»(PLC*Pl-P10) )/(2.-PLC> 


110 


112  P 


♦<l.-PLC)*(PlO-2.*Pin/(l.^PLC)^ 
12.*Pl7N*d,-.PLC)*tP12N-2,*P3Nn/(T,*PLC) 
2.*P18^{1,-PlCT*(P13N-2,*P5) )/ (I.^PlC) 


116 


120 


ii::?hE! 


PllNa(PLC*P16-2.*(PLC»P2-PllN) )/<2.-PLC> 
Pi2Na  (PLC*P17N-2,*  (PlC*P3N-P12N)  )  /  (2,->PL 


PLC) 


pi3Na(PLC*Pi8-2,*(PLC*P5-P13lQ)  )/(2.-PLC) 
P14Na(PLC*P19-2,*(PLC*P8-P14N) )/(2.-PLC) 

TU*APSO-UPSO 

TVaAPSQ-VPSQ 


82 


rinno 


TWaftPSQ-WPSQ 

CAPsC4^^1*TU*CAP2*TV*CAP3*UP*VP 

8cP*^CC^l*TU*PTP*(TV*£irP*CCP3*UP*VP>CCP4*VP*WP)  ♦CCP2*TW 
1  ♦CCP5*UP*WP) *HP 

C0P»CDPl*rv>CaP2*UP*VP 

CEP= (CEP1*TV*GTP*CEP2*UP*VP*CEP3*VP*WP) *HP 
CFP=  (CFPl*TU*CFP?*Ty*GTP«'UP*VP*(CFP3*GTP*CFP4) 

1  ♦WP* (CFP5*VP*CFP6*UP) ) *HP 

CJP=(CJP1*UP*VP*CJP2»TV)* (UP-CSA*CJP6*(WP-SNA) ) 

1  ♦ (CJP3*UP*VP*CJP4*VP*WP*CJP5*rv*HTP)*(WP»SNA) 

ACAP(K) »CCP*C3(K) 

BCAP(K)  aCAP*A2<  1  )*«X1^CBP«*82(  J)  ♦CCP«C2(K) 

CCAP(K>sCCP*Cl (K) 

OCAP (K) »CJP-CAP* ( A  1 ( I ) •P1>A2 ( I ) *«X2*P2*A3 ( I ) *P3N) -CBP* 

1  (B1  (J>*Pl6*B3(J)*PllN)-CDP*(Dl*(P15*P12N)*D2*(P10*P17:'i»  )- 

2  CEP*  <£1*(P19*P13N) ♦E2*(P18*P14N)) -CFP* (F I* (P7*P6N) ♦F2* (P4*P9N) ) 
IF  (K.EO.Kl)  UCAP(K)eOCAP(K)-ACAPlK)*P5 

IF  (K.E0.K2)  OCAP(K)«DCAP(K)-CCAP(K)*Pa 
IF  (I.LE.ITE.OR.J.LT.JROOT)  60  TO  700 

IF  (K.EQ.KW)  0CAP(K)*DCAP(K)*CCP*C1(K)*GAMMA(J)-C£P*E2* 

1  {GAMHA(J*1)-GAMMA(J-1) ) 

IF  (K.EQ.Kw^i;  UCAP(K)»0CAP(K>-CCP*C3<K)*GAMMA(J)-CEP*E2* 

I  (GAMMA(J>I)-GAMMA(J-ln 
60  TO  TOO 

DEFINE  POINTS  OF  HYPERBOLIC  (SUPERSONIC)  COMPUTATION  MOLECULE  AND 
COMPUTE  TRI-DIAGONAL  COEFFICIENTS 


200 


NSUPaNSUP*! 

P3N»P ( I-l ♦ J»K) 
P4*P(I*ltU»K-l) 

IF  (1,EvJ.1TE.AN0.K1.EQ.KW*1) 
P6NaP(I-l«J»K-l) 
P7*P(I*1»J»K*1) 

P9NaP ( I-l ♦ J*K>1 ) 
P15*P(I*1»J*1#K) 


P4«P4^6AMMA( J) 


(I»J*I»K4l) 
.,EU.l)  60 


P17NsP(I-l  _ 
Pl8aP(i»J»l.K* 
Pl9aP<I.J*ItK^ 

^F  (J.EU.l)  60  TO 


270 


P13NaP(If 

P14N*P(If 


270 


272 


276 


280 


h 


IF 

F 

0 

P 

P 

P 

P 


(J.EO.UK)  GO  TO  272 
(J.EU.JKP)  GO  TO  276 
TO  200 


03PSYM(I*1»K) 

IN«PSYM( I fK) 

2NsPSYM( I-l tK) 

3NaPSYM(I»K-l) 

Pl4NxPSYM( I ) 

GO  TO  280 

P15»(2.*P15^ (1.-PLC)*(P10-2.»PI) )/(l.*PLC) 
P17N*(2.*Pl7N* {l.-PLC)*(P12N-2i*P3N) )/(l,*PLC) 
P18»(2.*P18*(l,-PLC)*(P13N-2,*Pb) l/(l,*PLC) 
P^9*<2.JP19* ( 1,-PLC)*(P14N-2.*P8) ) /(l.*PLC) 

P12N*(PLC*P17N-2.*(PLC*P3N-P12N) ) / (2,-PLC) 
P13N«(PlC*P18-2.*(PLC*P5-P13N) )/ (2,-PLC) 

P14N» (PLC*P19-2.* (PLC*P8-P14N) )X (2.-PLC) 
CONTINUE 


P20»PT3(J»K) 
IF  (K.EQ.M) 

IF  (k.ne.k: 

IF  (K.EQ, 

IF  (K.NE, 


r.m 

.K2) 


P21*P5 

P,21=P(If  J*K-2) 
P22aP8 

P22aP( If JfK*2) 
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IF  (J.EQ.i)  P<;3*2.*Pl  1-P2 
IF  <J,£U.2)  P23*PSVMn»K) 
(J.GE.i)  - 


IF 


P23sPTUJ-2»K) 


CPP= 

COP* 

CPPs 

^CSP* 

CTP* 

CVP* 

I 

CWP» 

^UCP2 

UC02 

UCR2 

UCS2 

UCT2 

UCV2 

UCI»»2 

FEO 

CLP 

CMP* 

CNPa 


CAPl*UPSQ*CAP2*VPSQ-CAP3*UP*VP 

CBPl*VPSQ 

(CCP1*UPSU46TP*  <  VPSQ*GTP-CCPJ*UP*VP-CCP**VP*k<P)  ♦CCP2*ii<PSO 

-CCP5*UP*»IP)  *HP 

CDP1*VPSQ-CDP2*UP*VP 

(CEPl*VPSQ«GTP-CEP2*UP*VP-CtP3*VP*WP) •HP 
{CFPI*UPSQ»CFP2*VPSQ*6TP-UP*VP* (CFP3*ttTP^CFP4> 

(CFP5«VP*CFP6*UP)  >  •HP 

(CJP1*UP^VP-CJP2*VPSQ)*<UP-CSA-»CJP6*(»<P-SNA)  ) 

♦ (CJP3*UP^VP^CJP4*VP*WP-CJP5«VPSQ*HTP) • <WP«SNA) 

»QPSQ*(CAP1>CAP2)-CPP 

arOPSQ*CBPI-CQP 

»(3PSQ^{CCP1*GTP^GTP-»CCP2)^HP-C«P 

»QPS(i*CDPl-CSP 

aQPSQ*CEPl*6TP*HP-CTP 

sQPSO* (CFP1*CFP2«GTP)^HP-CVP  ^  ^  „ 

*QP5Q*  (-CJP2*  (UP-CSA*C  JP6»  <  ><P-SNA)  >  -C  JP5*HTP*  ( WP-SNA> )  -CWP 
FP*UP/ (CJ*OXI) 

FEQ^(UP/CsJ*VP^CB6> 

FEd^VP/BSPAN 

FEQ^ (UP^CU2*GTP*VP*WP*CW1) 


1 

2 


TQ*1.- 

ACAPf* 

aCAPT 


CCAPT 

OCAPT* 


1 

! 

4 

I 

1 


DCAPT* 


QPSU/APSQ 
i^jcp v#c 3  ( K  ) 

:yS5i:5  ’  ’ 

-EPSI«  (CLP*A3n)^DXI*CMP^B3tJ)*DETA» 

UCR2*C1 (K) 

-UCP2*<A1  <  n*<Pl-P2)  ♦A3{1)*P3N) 

-UCU2^(ai  (J)^(P16-P2)*83<'>>*PUN>  ^ 
-UCS2*(D1*(PI5^P12N)*D2^  PiO*Pl7N) 
-UCT2*(El*(Pl9^Pl3N>^E2*<Pla*Pl4N) > 

-UCV2* (Fl« (P7tP6NI ♦F2*(P4*P9N) )-UCW2 
ocapt-tq* 

(CPP*  <-A3  ( I )  •  (P2*P3N)  -A5i  I  >  *  <P3N-P20n 

♦  CQP«(-B3(J)*(P2*P11N)-8SU)^IP1IN-P23)  ) 

♦  4,*CSP*(01*P12N402*(P3N*PIINI  >  ♦C>^P> 
ocapt=ocapt*epsi^ 

I  {CLP*A3U)^DXI*(-P2-P3N*P3J  ,  ^ 

\  ♦CMP*B3(J)  *OETA*  (-P2-PllN^Pn)  > 

IF  (WP.EQ.O.)  GO  TO  320 
IF  iwP.LT.O.)  GO  TO  330 

ACAP(K)sACAPT*TQ*  t-CRP*(C3(K»*C5lK) ) (CTP*Ea*CVP*F2) ) 

♦£PSI*CNP»C3(«)*02ETA 

BCAPtK) sBCAPT*TQ*2,^(CPP*C3tK) ♦2.« (CTP*EI^CVP*Fl ) » 
•EPSi^CNP^C3(K)*02ETA 

|aO^APT-TQ»(CRP*(-C3(K)*P2*C5(K) •P2i) 

♦  4,*CTP*(El*Pl3N*E2*PnN> 

♦4.*CVP*(F1*P6N*F2*P3N) ) 
♦EPSI^CNP^C3(K)^02ETA*t-P2*P5> 

GO  TO  400 


AP<K] 
AP  (Kl 


320  ACAP(K 

^8CAP(K 

CCAP<K 

'oCAP(K 

1 

2 

GO  TO 

330  ACAP(K 
8CAP (K 
I 


|*ACAPT^TQ^(CRP*C3<K) ♦2»*<CTP^E2*CVP^F2) ) 
♦ePSi»CNP*HP/TOZETA 
>»BCAPf-TO*CRP*(Cl(K)*C3(«))  ^ 

) »CCAPT*TO^  <CRP*ci <K) ♦2.*ICIP*EI*CVP*FI) ) 
-tPSI*CNP*HP/rOZETA 
) aDCAPT-TQ* (2.*CTP« (EI*P13N^E2*P14N) 

♦  2,*CVP*(F1*P6N*F2»P9N>  » 

♦EPSI^HP* (-PS^PB) /T02ETA 
400 

>*ACAPT 

) sBCAPT*TQ*2,*  <CRP*Cl (K) *2. • (CTP*E2^CVP*F2) ) 
♦EPSl^CNP^Cl (K)^D2ETA 


84 


J 


noooo 


CCAP  (K)  sCCAPT^TO*  (-CRP*(C1(K»*CMK))  *4 ,  ♦  ( CTP*E  1  ♦CVP*F  1 )  ) 
1  -epsi*CNP*ci  (»<>*ozerA 

OCAP (K) cOCAPT-TQ* (CPP* (-C1 (K ) •PZ^C 7 ( K ) *P22) 

1  ♦4,*CfP*(F2*PUN«'ei*PUN) 

2  ♦4,*CVP* (F2«P9N*FI*P3N> > 

3  ♦£PSI*CNP*C1  <»<»«OZETA*  (P2-P8) 

400  CONTINUE 

IF  (K.EQ.M)  UCAP(K)»DCAP(K)-ACAP(K)*P5 
IF  (K.E0.N2)  UCAP(K)»OCAP(K)-CCAP(K)*P0 
700  CONTINUE 
RETURN 


750 


COMPUTE  TRI-OiAGONAL  COEFFICIENTS  IN  THE  TREFFTZ  PLANE 

ENTRY  TWIT 
GPaG{J) 

TC=TAU(U)*CH0H0(J) 

CV2»6P/0SPAN 

CW1=1,/TC 

00  900  K*K1,K2 
HP»M (K) 

GTPaCAPGTl ( JfK) ♦CAPGT2(IM,J) 

HTP»CAPhT1 (J»K) ♦CAPMTZ ( IM, J) 

P5«P(IM,J»K-1) 

P8«P( IM,J»K*1) 

IF  (J.EQ.l)  ^11sPSYM(IM»K) 

IF  (J.NE.l)  PllsPTl ( J-1«K) 

P16»P(IM.J*1»K) 

VP*CV2*(P16-PU  )  /TDETA-*GTP*MP#  (P8-Pb)/T0Z£TA 
WP»SNA*CWl*HP* (P8-P5) /TOZETA 
IF  (K.NE.KW.ANO.K.NE.KW^I )  GO  TO  750 
VPaVP-6TP*Hp*GAMMA ( J) /TOZETA 
WPa»fP-CWl*HP*GAMMA  ( J)  /TOZETA 

continue  ^ 

TVsAI2-VP*VP^„ 

C8Tatv*GPZ (BSPaN*8SPAN) 

rrTaHP* ( (AI2-WP*WP) / <  Tc*Tc) ♦OTP* ( TV*6TP-2t*VP*WP/TC) ) 
CtT*2,*GP*HP* (TV*GTP-VP«WP/TC) /BSPAN 
CJT*TC*  <WP-SNA) * (2.*VP*WP*CAPI0( J) -TV*HTP) 


800 


850 


P2»P( IMfJfK) 

PlSaP ( IM, J*l ) 
P19ap(lMf J^l 
IF  (J.EQ.l)  GO 
PllN3P(IM*J-l*K) 
P13N8P( IMf J-1 fK-1 ) 
Pl4NaP ( IMt J-1 .K*l ) 
GO  Tn  850 

PUNaPSYM(  IM.K) 
P13NaPSYM ( Im.K-1 ) 
PlANaPSYM ( IM.K* 1 ) 


flOO 


ACAP(K)aCCT*CJ(K) 

8CAP(K) aCdT *82 ( J)*RX1 *007*02 (K) 


1 


AP  (K) 


IF 


_JT-C8T* (81 ( J)*P16*82(J)*RX2*P2*B3(J)*P11N) 
-CET* (El* (P19*P13N) ♦EE* (P18^P14N) ) 

(K.EQ.Kl )  UCAP(K) =OCAP(K)-ACAP(K) *P5 


1 


(K.EQ.K2)  DCAP(K)=0CAP(K)-CCAP(K)*P8 
(J.LT.JROOT)  GO  TO  900 

-  - -  -'■T*CI(K)*6AMMA{J) 


900 


1  -CET*E2* (Gamma ( j* I ) -gamma (J-i) ) 
CONTINUE 
return 
END 


CT*C3(K)*6AMMA(J) 
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subroutine  INVERT(K1»K2) 

inverts  TRI-DIAGONAL  matrix  (VARSA*  matrix  iterative  anal..  P.195) 


DIMENSION  G6(^5) .WW(25) 


COMMON/ABCO  / 
COMMON/SOLVO/ 


ACAP(25J ♦8CAP<25) .CCAPJgS) »0CAP<25) 

DLl (49.19) .0L2(49.iv) .DSl (49) .052(49) .DUl (49.19) ♦ 
DU2(49.19) .0PMAX.UPMAXT.GAMMA(19) . IMAXI . I TERT . 
JMAXI .jMAXT.KMAXl.KMAXT.NSUP.P(49.19.a5) .PLU (19) . 
PNEW(25) .PN0SE.PSyM(49.25) .PTl (19.25) .PT2(19.25) ♦ 
PT3 (19.25) .PTEl (19) .PTE2( 19) .PWSL (49. 19) .PW8U(49*1 


CALCULATE  COEEFiCIENTb  OF  REDUCED  SYSTEM  OF  EQUATIONS 

(Kn  =CCAR  (Kl)  /RCAP  (K1  ) 

G6  (  K1  )  aUCAR  (Ki  )  /  BCAP  (KU 
KlPrKl*! 

DO  100  KrKlP,^2 
KZrK-1 

DENOMaBCAP  (K)-ACAP(K)*W(J(K7) 
rtW  (K  )  =CC/^^lt<)  /OENOM 

QG (K ) a (OCAP  <K ) -ACAP (K ) *G6 (KZ) ) /UENOM 

INVERT  reduced  MATRIX  BY  HACK  SUBSTITUTION 

PNEw (K2) =o6 (K2) 

KDaK2-K 1 

DO  200  Kal.KD 

KZ=K2-K 

PnER (KZ) aUG (KZ) -WW (KZ) •PNEW (KZ*1) 
return 

END 


SUBROUTINE  MAXI (Kl ♦K2. I , J) 


determines  location  and  absolute  value  of  the  maximum  change  IN 
POTENTIAL  between  THE  LAST  TWO  COMPUTED  ITERATES. 


COMMON/SOLVO/ 


DL 1(49. 19) .0L2(49.19) .DSl (49) .DS2(49) .001(49.19) ♦ 
DU2(49. 19) .OPMAX.OPMAXT.GAMMA(19) . IMAXI . I TERT. 
JMAXl .JMAXT.KMAXi.KMAXT.NSUP.P(49.19.25) .PLEl (19) . 
PNEW (25) .PN0SE.PSYM(49.25) »PT1 (19.25) .PT2( 19.25) ♦ 
PT3( 19.25) .PTEl (19) .PTE2(19) .PwBL(49.19) .PW9U(49*1' 


00  100  K=M,K2 
l)ELPaPNEW(K)-P(I,j,N) 
DPsARS (OELP) 

IF  (dp.lt.dpmax)  go  to 

DPMAXaDP 

IMAXIal 

JMAXIaJ 

KMAXfaK 

CONTINUE 

return 

ENO 
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nr>o  oorioor; 


SUHPOUTiNt  ARHAYS(NPRINT) 

OUTPUTS  SOLUTION  ARRAYS  EVERY  NPRINT  ITERATIONS.  IF  NPRINT.EQ.O. 
ONLY  THE  CIRCULATION  IS  WRITTEN. 


COMMON/CONST/  ALPHA»2MACW.0PLIM»|PSI;»wE.WS»BSPAN»AI2*CSA»RXI.RX2. 

•  SNA.TOETAtTOXI .TD/ET  A.OETAtOXI f OZETA* ILE * IM* INOSE » 

•  ITAIL  *ITE. JM.UHOOr.UTIP»KM»KW» jai« J82. J83fPLC8l. 

^  P  t- C  Q  2  •  PL.  C  9  3 

COMMON/SOLVO/  OLl  (49.19)  .DL2(49»19)  .OSl  *91? * 

*  0U2(49.19) «0PMAX.0PMAXT»GAMMA(I9) » IMAX I . I TERT. 

•  JMAXI » JMAXT.KMAxi.KMAXT  »NSUP .P ( 49 . 19 » 25) » PLE I ( 19 » . 

*  PnEW (25) * PNOSE ♦ PS YM( 49 .25) »PT1(19»25).PT2(19»25). 

*  PT3( 19.25) »PTE1 (19) ♦PTE2 (19) .PWBL (49. 19) .PWBU(49. 19) 


L 


C 


C 


c 


c 


WRITE  (6.161) 

WRITE  (6,171) 

write  (fa. 171) 

write  (6.181) 

WHITE  (6.161) 

WRITE  (6.131)  (GAMMA (J) ,J=l .JM) 

IE  (NPRINT. EO.O)  GO  TU  225 

write  (6.161) 

WRITE  (6.141) 

00  ISO  Jsi.JM 
write  (6.142)  J 
00  U5  Isl.IM 

14S  write  (6,146)  I ♦ <P< I . J.K) ,K=1 .KM) 
150  CONTINUE 

write  (6,161) 


105 

no 


WRI  ?E 
00  110 
WRITE 
00  105 
WRITE 


1 


(6.101) 
U=i  ♦  JM 
(6.102) 
1*1.  IM 
(6.106) 


CUNT INUE 
WRITE  (6.161) 


J 

I fPWBL ( I ,J) .OLl ( I .0) ,0L2 ( I ♦ J) .002(1,0) .UUl ( I, J) , 
PwPU( 1 » J) 


WRITE  (6.121) 

00  125  1=1. IM 

125  WHITE  (6,126)  I .nsi ( 1 ) ,052 ( 1 ) 

write  (6.161) 


220 

225 


WHITE 
00  220 
WHITE 

WHI  TE 

write 

WHI  Tf 
WHITE 

return 


(6,191) 
I  =  I  .  I M 
{6.14fa> 
(6,161) 
(b.lTl) 
(6.1/1) 
(6,161) 


I  f  (PSYM( I ,K)  ,K=1 ,KM) 


101  FORMAT 
1 

102  format 
I 

106  format 
121  format 
1 

126  format 
ill  format 
I 

141  format 
1 


(  1H-///1X40HDUMMY-PO1NT  POTENTIALS  ON  VERTICAL  LINES/ 
IX25H(PLANE-HY-PLANE  SPANWISE)) 

(///IXIOHPLANE  J  =  I3//3X1HI. 16X4RPWBL.12X3H0LI. 12X3H0L2. 

I  7X3HUU?f  12X3H0U1  .HXAHPWtiU/) 

(  I4,2<5X, 3€ 15.4) ) 

(1M-///1X3BH0UMMY-PO1NT  potentials  ALONG  SIDE  EDGE// 

3XIhI . ITXSHDSl .12X38052/) 

( l4,5X.2eiS,4) 

( 1M-///1X34HCIHCULATI0N  OlSTRiHUT ION,  6AMMA(J)// 
(4X.13E10,3)///) 

( ih-///ix26hp( I « J.K) ;  potential  array/ 

lX43H(PLANe''RY-PLANt  SPAN-.ISE.  I  OOWN.  K  ACROSS)) 
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r'.no  o  ooonn 


142  FOWMftT 
146  format 
161  format 
171  format 
1 
2 

181  format 

191  format 


(///IXIOHPLANE  J  =  13/) 
( 14/ (4A, 1  3E1 0. J) ) 
(//////) 


(  1H-///1X15RS0LUT10N  ARRAYS) 

(1H-///1X46HPSYM(I,K) !  IMAGE  POTENTIALS  AT  SYMMETRY  PLANE/ 
1X18H(I  OOWN,  K  ACROSS)/) 


SUBROUTINE  CPCOMP 

COMPUTES  the  WING/BODY  SURFACE  PRESSURE  DISTRIBUTION, 


OIMEN 

DIMEN 


SION  CPL<49» 
SION  ZXL(49. 


19) *CPU(49» 19) 
19) »2XU(49*19) 


COMMON/CAPS  / 
COMHON/CONST/ 

COMMON/COORD/ 

COMMON/SOLVO/ 

COMMON/SUHF  / 


^LPHA»|MAiHi6^LlM*i^s! ♦3E»Wb*BSPAN*Al2»CSA»«XltRX2» 
SNA,TO£TA,TOXI«TOZETA,DETA.OXI,DZETA#IlE*IM,IN6sE» 
ITAIL»ITE,JM,JR00T»JTIP*KM»XW» JB1*JB2*JB3*PLCBI* 
PLCB2  fPLCB3 

OLl (49«19) »0L2(49.l9) »DSl (49) *052 (49) .001(49* 19) ♦ 
DU2(49.l9) .OPMAX.OPMAXT. GAMMA (19) ♦ IMAX I ♦ I TERT » 
JMAXI.UMAXT.KMAXl»KMAXT.N5UP»P(49f 19*25) .PLEl (19) * 
PNEW(25) ♦PN0SE.PSyM(49,25) ♦PT1(19.25) »PT2(l9.25)  . 
PT3( 19,25) ,PTE1 (19) .PTE2(19) ,PWBL (49, 19) .PWBU (49. I' 


KB0DU(49.19) ,ZL (49. 19) .ZU (49* 19) 


♦0.?0/AI2)/1.20)**3.5-1,>/0.70 


CPSTAR=AI2* ( ( ( 1 
DO  100  1  =  1, IM 
DO  100  J=1.JM 
ZXU ( I , J) =0. 

ZXL ( I ,0) =0, 

CPU(I,U)=0.0 
100  CPL(I,U)=0.0 

DK=DETA 

IMX=IM-1 

00  200  MSURF=1,2 

DO  200  IsINOSE.IMX 

JEsJPOD ( I ) -1 

00  200  JsI.JE 

IF  (MSURF.EU.2)  GO  TO  110 

t<P  =  -l 

OL=-nZETA 

UEL=-0ELL (I , J) 

I2=IZL (J) , 

KS=KBOOL ( I . J) 

^  ^  ^  9  ) 

CAPGTB=  (OL-DEL)  •  (CAP{-.ri  ( J,KS) -CAPoT  1  (J.KS*!)  )/DL4CAPGTl  (J,KS*1) 
CAPGTB=CAPGTB*CAPGT2( I ,J) 

GO  TO  120 
110  KP=*1 
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c 


c 


c 


■>L=07ETA 
i3bL=LlELU  (  I  .  J) 

KS=NBOOU ( i • J) 

H 5H WBU ( i * J ) 

CAPCjTW  =  (OL-DEL  )  *  (CaEB  [  1  ( J.KS)  -CA^'OTl  IJ,KS-1  )  )  /DL^CAPGT  1  ( J»Kb-l ) 
CAEbTB  =  CAPijTrt*CAPGT?  (  I  ♦  J) 
l,an  IF  (J.GE.JHUOT)  GO  Tn  140 
DHs-OXi 
IP  =  -1 
GO  TO  IHO 

140  IF  (I.GE.IJ)  GO  TO  160 
UHs-OXl 
IP  =  -1 

60  TO  180 
160  DH=*DXl 
IP=>1 


180  Al=-1 .b/Db 
A2=2.0/Drt 
A3  =  “n,6/nri 
81=-1 .B/OG 
B2  =  /?.0/OK 
B3  =  -(l  .b/Uix 

Cl=-  (I3L*?.*1)El  )  /  (PEL*  (OL*OFL3  ) 

CB=  <OL*D£L) / (uEL*Ol  ) 
C3=-nEL/(OL*<OL*OEL) ) 

01  =  0,5*  (DL*  DEL)  •  (  2  .  *nL*DEL )  /  ( DL*t)L  ) 
02  =  -nEL*  (^l.*nL*0EL)  /  (uL*OL) 
Q3=0.5*0El* (UL*nEL) / (UL*OL) 


FP  =  E  (  I ) 

6P=6( J) 

CJ=CH0HD (0) 
CHG=CAP6( I . J) 
CUlsFP/CJ 
CU2=CAPGy ( I . J) 
CVlsCBO*FP 
Cv2=GP/8SBAN 
CWl=l ,/ (TAU ( J) -CJ) 


181 


1  18? 


18? 

183 


P2=PWHL ( I . J) 
P2  =  P»NrfO  ( I » J) 


I3=l4lP 
I4  =  I3MP 
J5=J4 1 
J6=J*? 

KS1=KS*KP 
Kb2=KSl*KP 
IE  (MSURF.EU.l) 

IF  (MSURE.EQ,2) 

Pl=P<I»Jf^S) 

P?  =  P ( I  ♦  JtGSl ) 

P3aD]*PlIi*Jt^5)*0?*P(I3»J»‘<Sl)*03*P(I3,JfKS2) 

P4si31*P(14,J,I^5)*D2*P<IA,J,KSI)  *03*8  114  t  J*KS2) 

Pb  =  l3l*P  (  I  I  J5.^S)  ♦OE^P  (  I .  JSfKSP  ♦U3*P  ( I  ,  J5*KS2) 

P6  =  i31*P  (  I  t  Jb.tXS)  ♦f3?*P  (  I  ♦  JbtKSl  )  ♦U3*P  (  1  ,  J6»KS2) 

IF  (  J.EO. Jfll-1 .OP.J.FO, JBl )  GO  lu  181 

IF  (  J.EQ, JH2-1 .0  ?, J.EO. JB2)  GO  TO  1182 

IF  ( J,Eu. JH3-1 .OP. J.EO. JB3)  GO  10  182 

GO  TO  118b 

Jt<  =  J81 

PLC=PLCbl 

GO  T(i  183 

JK=JH2 

plC=PLC82 

GO  TO  183 

JK=JB3 

plc=plcb  i 

COnT INUE 

PA=(?.*P(i«J*^*l*KS)*(l  ,-PLC)  *  (P  <  I » JIX-1  *KS) 

1  -2.*P < I «0K ,Kb ) ) ) / ( 1 . ♦PLC)  ,  , 

PH= (?.*P ( i t JK* 1 .KSl ) ♦ ( 1 .-PLC) *  IP ( I » JK-l.KSl ) 
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c 


1  -2.*P(1.JK,KS1) ) )/(l.*PLC) 

PC=  (2.*P  ( I»  JK^l  ,KS2)  ♦  ( 1  .-PLC>«<Pn  f  JK-1  «KS2) 

1  -2.*P<I»JK,KS?) ) )/(l.^PLC) 

IF  (J.EQ.JK)  GO  TO  184 
P6=0l*PA>02*Pb^D3*PC 
GO  TO  1185 

184  P5=U1*PA*U2*PH*D3*PC 

PAA=P ( I » JK-1 ♦KS) +3.* (PA-P(I ♦ JK*KS) ) 
P8B=P(l»JK-l»KSl) ♦3.*(P8-P(I*JK*KSl) ) 
PCC=P(I.JK-1*KS2)  ♦G.^CPC-Pdf  JK*K!j^)  ) 

P6  =  l}l*PAA*02*P8BtD3*PCC 
1188  0PDS=A1*PZ*A2*P3*A3*P4 
DPD2=CI*PZ^C2*P1>C3»P2 
u=csA^cui*opns>cu2*HWri*opnz 

V=CV1*DP0S>CV2* (B1*P7*R2*P5*Q3*H6) ♦CAP6TB*hW0*DPDZ 
y<sSNA*C41*HVi8«0P0Z 
IF  (J.EQ.l)  V=0.0 
OX=l ,-u*u-v*v-w*w 

CALCsAIZ* I ( 1 .♦0.20*QX/A12)**3.5-1. ) /0.70 
IF  (MSUHF.EU.Z)  GO  TO  185 
ZXL  ( I  ♦  J)  =  (W-V*OZOYL  ( I  ♦  J)  )  /IJ 
CPL<NJ)=(-ALC 
GO  TO  200 

185  ZXU(I»J)=(W-V*DZDYU(I*J) )/U 
CPU(t»J)=CALC 

200  CONTINUE 


1 

310 


JROOTX=JMUOT-l 
WRITE  <6»J01)  CPSTAR 
00  310  J=1»JR00TX 
YY=2,*Y (J) /BSPAN 
WRITF  (b»302)  J»YY 
00  310  IslNOSE.IMX 
TEMPaXCAP(I) ♦O.B 
OZOXlaUZOXUdfJ) 

0Z0X2=0Z0XL ( I tJ) 

WRITE  (6f303>  T£MP,CPUd»J)  »2XUd.J)  »L)ZDXl*CPLd*J)  » 
ZXLdtJ)  «0Z0X2 

CONTINUE 

00  320  JsJR00T«JTIP 
YY=2.*Y (J) /BSPAN 
WRITE  (b»302)  JtYY 
00  320  IsiLEdTE 
TEMPsXCAPd)  ♦O.B 
OZDXlsUZOXUdfJ) 

L)ZDX2=D20XLd»J) 

WRITE  (6t303)  TEWP » CPU d ♦ J) t ZXU <1  * J) * DZDX 1 fCPL d * J> f 


320  CONTINUE 


ZXLdfJ)  f0Z0X2 


return 

C 

301  FORMAT 
1 

30?  format 
1 
2 

303  FORMAT 


dH-//////lX30HSURFACE  PRESSURE  DISTRIBUTIONS/// 
1X8HCPSTAR  =  E13.5///) 

dri-fWHSPAN  STATION  J  «  I4*6X6HZY/e  =  F9,4// 
3X9m(X-xlE» /Cf 12X3HCPU*8XTHSL0PE  U» 10X5HDZDXU» 12X3HCPL* 
8X7HSLOPE  LflOXSHOZOXL/) 

(Fi2,4»2(F15.4,2Fl5.6)  ) 


END 
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